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1.  Project  Summary 


In  this  project  we  investigated  the  behavior  of  early  body- wave  coda  observed  at  intermediate 
distances  from  events  in  Central  Asia.  At  these  distances  (~14°-  29°  degrees)  the  waveforms 
contain  triplicated  arrivals  from  upper-mantle  discontinuities  that  could  potentially  improve 
seismic  monitoring  capabilities.  However,  intermediate-distance  observations  are  typically 
under-utilized  in  location  and  magnitude  estimation,  because  heterogeneity  along  the  path  and 
complicated  phase  interactions  make  the  body- wave  arrival  suite  difficult  to  interpret. 

The  phase  complexities  that  occur  at  distances  between  13  -  30°  produce  uncertainties  and 
errors  in  the  phase  picks  in  seismic  bulletins.  We  demonstrate  this  phenomenon  in  Figure  1,  in 
which  we  plotted  Pn/P  travel-time  residuals  as  a  function  of  epicentral  distance.  To  generate  this 
figure,  we  retrieved  International  Seismic  Centre  (ISC)  bulletins  for  the  years  1999-2001  from  the 
Flinn-Engdahl  Seismic  Regions  26-30,  47  and  48  (see  http://neic.usgs.gov/neis/epic/fer.htmlT 
which  include  most  of  Asia  and  parts  of  the  Middle  East.  We  included  all  events  with  ISC  depths 
less  than  40  km  that  were  located  by  more  than  25  stations.  In  Figure  la,  we  plotted  the  residuals 
(gray  dots)  with  respect  to  the  IASPEI91  (Kennett  and  Engdahl,  1991)  model  for  associated 
first-arriving  P- wave  picks  between  10  -  30°.  Then  we  calculated  the  density  (number)  of  residual 
picks  in  boxes  of  0.5°  distance  by  0.5  seconds  of  residual,  and  superposed  a  smoothed  version  of 
the  resulting  image  over  the  individual  picks.  We  also  plotted  two  dashed  black  lines  in  Figure  la 
that  show  the  travel-time  differences  between  branches  C’B’-  A’B ’  and  CB  -  AB,  with  the 
corresponding  ray  diagrams  and  travel-time  branches  shown  in  Figure  lb. 

a)  b) 


Epicentral  Distance  (degrees) 


Epicentral  Distance  (deg) 


Figure  2.  a)  P/Pn  phase  residuals  (with  respect  to  the  IASPEI91  model)  from  ISC 
bulletins  in  Asia  (1999-2001)  as  a  function  of  epicentral  distance;  the  residual  density  as  a 
function  of  distance  and  reduced  travel  time  is  shown  as  a  superposed  color  image.  Heavy 
black  dashed  lines  show  the  travel-time  difference  between  secondary  and  primary 
upper-mantle  phase  branches  of  the  iasp91  travel-time  curves,  shown  in  b). 


1 


Figure  la  clearly  shows  a  large  cloud  of  travel-time  residuals  that  are  associated  with  the 
410-km  triplication  (occurring  at  14  -  18°  distance).  There  is  also  considerable  structure  to  the 
residuals  in  the  regional  to  far-regional  range,  revealed  by  the  negative  residual  bias  at  distances 
between  10  -  17°.  While  this  particular  data  set  does  not  reveal  an  increase  in  phase  residuals 
associated  with  the  660-km  discontinuity  (i.e.,  the  CB  travel-time  branch),  we  have  observed  such 
an  effect  for  other  time  ranges  and  in  other  bulletins  (e.g.,  National  Earthquake  Information  Center 
[NEIC]).  This  type  of  residual  behavior  as  function  of  epicentral  distance  is  also  present  in  the 
phase  picks  used  to  locate  events  at  the  ISC,  but  other  agencies  restrict  their  location  travel-time 
picks  more  aggressively,  making  the  phenomenon  shown  in  Figure  la  less  apparent. 

The  large  arrival-time  residual  patterns  displayed  in  Figure  1  provided  the  original  motivation 
for  our  research  project.  The  pattern  of  the  residuals  suggested  that  multi-pathing  effects  from  the 
410-km  and  660-km  discontinuities  causes  phase  arrivals  to  be  systematically  mispicked  and/or 
mislabeled.  Triplication  effects  from  multi-pathing  produce  amplitude  and  arrival  phenomena  that 
result  in  uncertainties  and  errors  in  the  phase  picks  in  seismic  bulletins. 

To  address  these  issues,  we  developed  a  set  of  characterization  techniques  for  phases  observed 
at  intermediate  distances,  combined  with  studies  of  Earth  structure  along  paths  sampling 
upper-mantle  transition  zone  structure  in  Central  Asia.  We  have  applied  our  techniques  to  data 
from  two  arrays  in  Kazakhstan:  MKAR  (Makanchi)  and  KKAR  (Karatau),  both  of  which  record 
far-regional  events  from  most  of  south-central  Asia  (see  Figure  2).  To  provide  a  test  bed  for  our 
analysis  techniques,  we  collected  waveform  data  for  moderate-sized  events  observed  at  the 
MKAR  and  KKAR  arrays.  We  populated  the  database  with  events  that  were  well  located 
teleseismically  and  had  Global  Centroid  Moment  Tensor  (CMT)  solutions  (http://globalcmt.org) 
associated  with  them.  In  total  we  collected  waveforms  on  600  events  that  appear  in  the  EHB 
bulletin  (Engdahl  et  al.,  1998)  for  the  years  2002  through  2006.  Figure  2  shows  the  EHB  locations 
for  these  events,  which  have  epicenter  distances  ranging  between  14°  and  29.7°  from  at  least  one  of 
the  arrays.  Most  of  the  events  are  located  in  the  Middle  East  and  Southern  Asia,  extending  through 
Mongolia  to  southern  Siberia.  There  are  no  events  to  the  north-northwest  of  the  arrays  in  the 
database,  due  to  lack  of  seismicity  in  these  regions  at  the  intermediate  distance  range.  The  database 
contains  296  events  recorded  at  MKAR  and  304  events  recorded  at  KKAR,  with  an  additional  118 
events  that  are  in  the  intermediate  distance  range  at  both  KKAR  and  MKAR.  Magnitudes  vary 
between  mb  =  4.0  -  6.8,  with  most  events  below  an  mb  of  5.8.  Events  depths  are  mostly  constrained 
to  the  crust,  but  there  are  some  events  with  hypocentral  depths  up  to  150  km  deep.  There  are  only 
38  events  in  the  database  that  have  depths  greater  than  38  km. 
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Figure  3.  Map  of  south-central  Asia  showing  the  location  of  the  Makanchi  (MKAR)  and 
Karatau  (KKAR)  arrays  in  Kazakhstan,  as  well  as  the  location  of  ~500  earthquakes 
(circles)  used  in  this  study.  Colored  bands  mark  the  14  -  29°  distance  range  from  each 
array. 

This  report  is  organized  in  the  following  fashion:  we  first  provide  an  in-depth  description  of  the 
seismicity  and  tectonics  that  are  observed  in  eight  distinct  zones  at  far-regional  distances  from  the 
MKAR  and  KKAR  arrays  in  Kazakhstan.  Next,  we  describe  the  difficulties  of  far-regional  phase 
analysis  observed  on  the  small  apertures  of  the  MKAR  and  KKAR  arrays.  In  the  fourth  section  of 
the  report  we  show  the  results  from  our  application  of  improved  array  signal-processing 
algorithms  that  are  specifically  designed  for  small-aperture  arrays. 

In  the  remainder  of  the  report,  we  describe  the  methods  we  developed  to  better  understand  the 
path-specific  structure  that  affects  seismic  body  waves  observed  at  far-regional  distances  from 
MKAR  and  KKAR.  For  example,  in  the  fifth  section  of  the  report  we  used  our  array-processing 
results  (slowness  and  phase-arrival  times)  to  empirically  derive  regionalized  velocity-depth 
profiles  that  more  accurately  predict  the  far-regional  phase  succession.  These  1-D  velocity-depth 
models  are  based  on  x-p  transformation  of  both  array  measurements  and  waveforms  to  form 
generalized  representations  of  the  path-specific  earth  structure. 

In  the  sixth  section  of  the  report,  we  describe  a  waveform  clustering  algorithm  that  we  devised 
to  explain  phase  behaviors  that  are  not  predicted  by  the  regionalized  velocity-depth  profiles.  The 
clustering  algorithm  first  sorts  a  set  of  array  beams  to  derive  'wavefield  templates';  i.e.,  grouped 
observations  with  similar  phase  characteristics.  Then  the  wavefield  templates  are  further  analyzed 
in  a  non-linear  fashion  by  comparing  them  with  synthetic  seismograms  that  are  derived  using  the 
regionalized  1-D  velocity-depth  models. 
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In  the  final  section,  we  describe  two  methods  we  used  to  examine  near-array  earth  structure 
and  its  effects  on  array  measurements.  This  involved  the  analysis  of  back-azimuth  residuals  from 
polarity  studies,  receiver  functions  to  image  below  the  arrays,  and  teleseismic  PcP  arrivals.  In  an 
appendix  we  demonstrate  how  small-aperture  arrays  such  as  MKAR  and  KKAR  might  be 
improved  for  far-regional  phase  analysis  studies  through  the  judicious  use  of  sub-arrays  or  the 
addition  of  some  carefully  placed  additional  array  elements. 

This  work  has  resulted  in  a  methodology  that  improves  the  phase  characterization  of 
far-regional  earthquakes  observed  on  regional  small-aperture  array.  Our  research  has  also  yielded 
insight  into  body-wave  phases  that  are  regularly  observed  on  the  MKAR  and  KKAR  arrays, 
including  information  on  expected  wave  propagation  behavior  and  the  regional  nature  of  the 
upper-mantle  discontinuities. 


2.  Far-Regional  Seismicity  Observed  at  the  KKAR  and  MKAR  Seismic  Arrays 

In  the  initial  phase  of  the  project  we  gathered  relevant  background  information  on  the 
seismicity  and  tectonics  at  far-regional  distances  from  the  MKAR  and  KKAR  arrays.  Earthquakes 
occurring  in  these  regions  exhibit  considerable  complexity  in  the  recorded  seismograms,  due  to 
differences  in  source  mechanism,  earthquake  depth,  near-source  structure,  and  upper-mantle 
discontinuity  structure.  Our  specific  goal  was  to  regionalize  the  earthquake  depths  and  source 
mechanisms,  which  gain  insight  into  the  phase  behavior  observed  in  our  database  of  MKAR  and 
KKAR  waveforms.  To  facilitate  this  analysis  we  divided  the  seismic  regions  into  separate  zones  of 
roughly  similar  tectonic  regimes  (shown  in  Figure  3),  each  of  which  we  discuss  in  the  following 
subsections. 


2.1  Zone  1:  Iran  and  Caspian 

We  derived  our  information  on  this  zone  of  seismicity  from  several  recent  studies,  including 
Tatar  et  al.  (2004),  Talebian  and  Jackson  (2004),  and  Engdahl  et  al.  (2006).  Most  of  the  seismicity 
in  Iran  results  from  the  collision  between  Eurasian  and  Arabian  plates,  forming  a  1000  km- wide 
zone  of  compression,  which  includes  continent-continent  collision  and  subduction  (Figure  3). 
While  the  Zagros  fold-and-thrust  belt  of  SE  Iran  is  the  most  seismically  active  region  of  Iran,  it 
falls  outside  the  14  -  18°  distance  range  we  focused  on  during  the  project;  thus,  it  was  not  included 
in  our  study.  Following  previous  researchers,  we  divided  the  region  labeled  as  Zone  1  in  Figure  3 
into  the  tectonic  provinces  shown  in  Figure  4.  These  provinces  include  1)  the  Caspian  Sea  region 
including  the  Alborz  and  the  Talesh  range;  2)  the  Oman  line  and  the  Makran  range;  and  3)  eastern 
Iran,  which  we  discuss  below. 

The  band  of  earthquakes  that  crosses  the  central  Caspian  Sea  (Aphsheron-Balkhan  sill)  is 
thought  to  represent  subduction  of  the  South  Caspian  basin  beneath  the  central  Caspian  (Jackson  et 
al.,  2002).  The  south  Caspian  basin  lacks  seismicity  and  likely  behaves  as  a  rigid  block. 
Earthquakes  across  the  central  Caspian  occur  to  depths  of  80  km,  and  there  is  little  evidence  that 
earthquakes  shallower  than  30  km  occur  on  the  offshore  portion.  Focal  mechanisms  (shown  in 
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Figure  5)  indicate  normal  faulting  with  the  T-axis  oriented  N  to  NNE.  Some  strike-slip 
mechanisms  have  also  been  reported  (Jackson  et  al.,  2002).  On  shore,  shallower  earthquakes  occur 
(30  -  40  km),  and  focal  mechanisms  show  thrust  faulting.  It  is  thought  that  the  normal  faulting 
earthquakes  do  not  represent  interaction  of  the  Eurasian  and  Arabian  plates,  but  rather  plate 
bending  and  slab  elongation  from  subduction. 


20° 


Figure  4:  Far-regional  seismicity  (14  -  18°)  from  the  KKAR  (red  circles)  and  MKAR  (blue 
circles)  arrays.  Other  seismicity  is  shown  in  yellow.  The  seismic  zones  discussed  in  the 
text  are  outlined  and  numbered.  Epicenter  locations  are  from  the  EHB  catalog  of 
Engdahl  et  al.  (1998).  Zone  1:  Iran  and  Caspian,  Zone  2:  southern  Pakistan,  Zone  3: 
northern  Pakistan  and  Afghanistan,  Zone  4:  Tibet,  Zone  5:  southern  Tarim  Basin,  Zone 
6:  western  Mongolia,  Zone  7:  central  Mongolia,  Zone  8:  Southern  Siberia. 

Along  the  south  Caspian  Basin,  earthquakes  in  the  Alborz  mountains  occur  shallower  than  15 
km,  although  events  in  the  EHB  catalog  show  depths  of  ~30  km.  Most  seismicity  here  exhibits 
left-lateral,  strike-slip  faulting  or  thrust  faulting,  with  trends  parallel  to  the  regional  strike.  A  few 
thrust  mechanism  strike  perpendicular  to  the  regional  trend  and  may  be  associated  with 
termination  of  strike-slip  faults  (Jackson  et  al.  2002).  Overall,  the  mixture  of  strike-slip  and  thrust 
faulting  likely  accommodates  the  oblique  left-lateral  compression  between  central  Iran  and  the 
southern  Caspian. 

The  west  coast  of  the  Caspian  includes  the  N-S  trending  Talesh  range,  a  continuation  of  the 
Alborz.  Seismicity  here  is  deeper  than  along  the  Alborz,  ranging  between  15-26  km  depth.  Focal 
mechanisms  for  these  earthquakes  (see  Figure  5)  indicate  very  shallow  angle  thrust  faults  (near 
horizontal),  with  slip  vectors  oriented  to  the  west. 
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Figure  5:  Seismicity  in  Iran  and  the  Caspian  Sea  region  (from  Engdahl  et  al.,  1998).  Red 
circles  are  earthquakes  observed  at  14  -  18°  from  KKAR  array.  Gold  circles  show  the 
seismicity  of  the  Zagros  and  other  earthquakes  of  the  region.  Localities  mentioned  in  the 
text  are  labeled  as  follows:  Zagros  (Z),  Oman  line  (OL),  Makran  (M),  Lut  block  (L), 
Central  Iran  block  (C),  Alborz  (A),  and  Talesh  (T). 

The  Oman  line  marks  the  transition  between  the  Zagros  continent-continent  collision  and  the 
subduction  of  the  Arabian  plate  beneath  the  Makran  coast.  A  change  in  earthquake  depth  is 
observed  from  the  shallow  upper  crustal  seismicity  of  the  Zagros  (<  20  km)  to  the  subcrustal 
seismicity  in  the  Makran  region.  Earthquakes  deepen  northward,  ranging  from  8-12  km  near  the 
coast  to  a  depth  of  28  km  just  50  km  to  the  north.  Most  focal  mechanisms  show  thrust  faulting, 
although  some  strike-slip  mechanisms  appear  in  the  Harvard  CMT  catalog  (Figure  5).  In  the 
Makran  region,  seismicity  is  seen  throughout  the  crust  and  extends  into  the  upper  mantle  to  depths 
greater  than  50  km  (Maggi  et  al.,  2000).  This  deeper  seismicity  likely  occurs  on  the  down-dip 
portion  of  the  subducting  Arabian  plate.  Focal  mechanisms  for  the  deeper  seismicity  show  normal 
faulting  with  T-axis  oriented  to  the  north  and  some  N-S  strike-slip  faulting. 

In  eastern  Iran  most  seismicity  occurs  along  the  boundaries  of  the  Lut  and  central  Iranian 
continental  blocks,  which  are  relatively  stable  and  aseismic.  This  region  has  historically  produced 
large  and  devastating  earthquakes  (e.g.,  Tatar  et  al.,  2004).  Earthquakes  occur  mostly  in  the  upper 
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crust  above  20  km  depth.  The  EHB  catalog  shows  a  mean  depth  of  12±5  km  (Engdahl  et  al.,  2006). 
On  the  western  edge  of  the  Lut  block,  the  focal  mechanisms  shown  in  Figure  5  indicate 
predominantly  NNE-SSE  strike-slip  faulting.  These  faults  likely  accommodate  the  right-lateral 
shear  between  central  Iran  and  Afghanistan  to  the  west.  On  the  eastern  boundary  of  the  Lut  block, 
a  NNW-striking  thrust  fault  appear  to  be  the  dominant  mechanism;  however,  strike-slip  faults  are 
also  apparent  in  the  Harvard  CMT  catalog. 

We  also  compiled  information  regarding  the  crustal  structure  across  the  region.  Most  of  this 
information  is  extracted  from  the  CRUST2.0  model  (Bassin  et  al.,  2000).  Some  information  is  also 
gathered  from  receiver  function  studies  (e.g.,  Mangino  and  Priestley,  1998).  The  various  regions 
show  large  variation  in  crustal  thickness.  Across  the  southern  Caspian  basin  the  crust  ranges  from 
30  to  over  45  km  thick,  with  the  thickest  crust  onshore  Across  the  Zagros  to  central  Iran,  the  crustal 
thickness  increases  from  45  km  thick  to  over  70  km  thick,  showing  an  abrupt  increase  east  of  the 
main  Zagros  thrust  (Paul  et  al.,  2006).  Profdes  extracted  from  CRUST2.0  indicate  a  variation  of  3 1 
-  39  km  in  crustal  thickness. 
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Figure  6:  Focal  mechanism  for  Iran  and  the  Caspian  Sea  region  discussed  in  the  text. 
Compressional  quadrants  are  darkened.  The  best  double-couple  Harvard  CMT  solutions 
are  shown  in  black.  Blue  mechanisms  are  from  Jackson  et  al.  (2002)  for  the  Caspian 
region.  Mechanisms  for  the  Zagros  are  not  shown. 
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2.2  Zone  2:  Southern  Pakistan 


Seismicity  in  southern  Pakistan  (Figure  6,  left)  is  at  far-regional  distances  (14  -  18°)  from  the 
KKAR  array.  Most  information  presented  here  is  summarized  from  Quittmeyer  and  Jacob  (1979) 
and  Ambraseys  and  Bilham  (2003).  The  majority  of  earthquakes  in  southern  Pakistan  result  from 
the  convergence  between  the  Eurasian  and  India  plate.  In  Pakistan,  high  seismicity  is  observed 
along  the  Kirthar  and  Sulaiman  fold-and-thrust  belts,  which  form  the  western  boundary  of  the 
collision.  Seismicity  also  occurs  along  the  Makran  mountain  range,  along  the  south  coast. 
Earthquakes  in  the  Sulaiman  range  fall  outside  of  the  far-regional  distance  focus  of  this  study  and 
were  not  part  of  our  analysis. 


Figure  7.  Left:  Reference  map  of  Pakistan  and  Afghanistan  with  places  discussed  in  the 
text  labeled.  Right:  Map  showing  regional  seismicity,  red  and  blue  circles  are  at 
far-regional  distances  from  KKAR  and  MKAR,  respectively.  Gold  circles  show  other 
regional  seismicity.  Focal  mechanism  are  the  best-double  couple  solution  from  the 
Harvard  CMT  catalog. 

The  Kirthar  range  is  bounded  on  the  west  by  the  Chaman  fault,  a  1000-km  long  left- lateral 
fault  system  that  likely  marks  the  boundary  between  the  Indian  plate  to  the  east  and  the  Eurasian 
plate  (Ambraseys  and  Bilham,  2003).  From  the  Harvard  CMT  catalog,  focal  mechanisms  show 
strike-slip,  thrust  and  normal  faulting  within  the  Kirthar  region  (Figure  6,  right).  The  normal 
faulting  events  occur  at  depths  greater  than  50  km  and  may  be  related  to  the  subduction  of  the 
Arabian  plate  beneath  the  Makran  coast.  In  general,  earthquakes  with  well-determined  depths 
occur  above  30  km,  at  an  average  depth  of  18  km,  similar  to  other  continental  regions.  The 
shallower  earthquakes  show  left-lateral  mechanisms  that  strike  NNE-SSW  and  thrust  faulting  with 
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the  T-axis  oriented  NW,  consistent  with  strain  partitioning  due  to  oblique  convergence  between 
the  Indian  and  Eurasian  plates. 

We  have  not  been  able  to  find  peer-reviewed  studies  of  the  local  crustal  structure  (e.g.,  receiver 
functions).  However,  the  CRUST2.0  model  (Bassin  et  al.,  2000)  shows  a  trend  of  crustal 
thickening  from  the  Makran  coast  north  into  the  Kirthar  range,  where  the  thickness  of  the  crust  is 
on  the  order  of  45  km. 

2.3  Zone  3:  Northern  Pakistan  and  Afghanistan 

Earthquakes  in  this  zone  are  at  far-regional  distances  from  the  MKAR  array  (Figure  7).  The 
region  includes  the  western  Himalayan  syntaxis  beginning  in  northern  Pakistan  and  extends 
through  the  Pamir  Hindu  Kush  region  of  Afghanistan.  Seismicity  in  Tajikistan  and  Uzbekistan  are 
also  included  in  this  zone.  The  main  structural  features  trend  northwest-southeast,  and  most 
seismicity  is  associated  with  Indian  and  Eurasian  collisional  plate  boundary. 
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Figure  8.  Seismicity  in  northern  Pakistan,  Pamir  Hindu  Kush  and  Uzbekistan  (from 
Engdahl  et  al.,  1998).  Blue  circles  are  earthquakes  observed  at  far-regional  distances 
from  the  MKAR  array. 

In  northern  Pakistan  the  main  seismic  zone  trends  northwest  following  the  trend  of  the 
Himalayas.  Source  mechanisms  show  mostly  thrust  fault  motion  with  NW-SE  oriented  fault 
planes.  Some  strike-slip  mechanisms  are  also  observed.  The  earthquakes  tend  to  have  shallow 
focus  (confined  to  the  mid-and-upper  crust)  and  most  are  associated  with  the  northward  dipping 
fault  system  that  separates  the  India  plate  from  the  Eurasian  plate.  Receiver  function  studies 
(Zugibe  et  al.,  2008)  indicate  that  the  crust  in  this  region  is  on  the  order  of  60  km  thick  along  the 
Himalayan  front. 
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The  Pamir  Hindu  Kush  seismic  zone  is  an  active  region  of  intermediate  depth  earthquakes  and 
has  seismicity  that  extends  to  300  km  depth  (Searle  et  al.,  2001).  The  mountain  ranges  are 
relatively  old  and  are  thought  to  have  been  reactivated  during  the  Indian-Eurasian  collision.  The 
convergence  of  the  Pamir  range  with  Eurasia  also  causes  shallow  focus  earthquakes  (<  45  km) 
along  the  thrust  fault  system  between  the  two  plates.  The  Hindu  Kush  seismic  zone  is  bounded  on 
the  east  by  the  Chaman  left-lateral  fault  system  (Figure  8)  and  on  the  west  by  the  right-lateral 
Karakorum  fault  (Searle  et  al.,  2001). 

The  pattern  of  seismicity  in  the  Pamir  Hindu  Kush  seismic  zone  is  quite  complicated  at  depth 
(Billington  et  al.,  1977;  Roecker,  1982).  A  large  scale  relocation  study  (6000  events;  Pegler  and 
Das,  1998)  indicates  that  this  seismic  zone  has  an  S-shaped  pattern.  At  depth,  the  earthquakes  map 
a  700-km  long  zone  that  is  30  km  wide,  which  potentially  represents  thin  subducted  continental 
crust.  It  has  been  hypothesized  that  the  northward  dipping  seismic  zone  occurs  on  Indian 
lithosphere  subducting  beneath  the  Hindu  Kush,  and  the  southward  dipping  seismicity  occurs  on 
Asian  lithosphere  underthrusting  the  Pamir  range  (Burtman  and  Molnar,  1993). 

North  of  the  Pamir  Hindu  Kush,  seismicity  in  the  Tajikistan/Uzbekistan  region  occurs  in  an 
intraplate  setting,  where  earthquakes  are  occurring  in  Asian  lithosphere.  Strain  in  the  region  is 
associated  with  the  Indian  and  Eurasian  collision  and  many  thrust  faults  exists,  although  strike- 
slip  also  occur.  Most  earthquakes  show  thrust  mechanisms  with  various  fault  plane  orientations.  A 
study  of  the  Gazli  region  of  Uzbekistan,  where  magnitude  7  earthquakes  have  occurred,  indicates 
that  most  earthquakes  occur  above  25  km,  likely  marking  the  limit  of  the  brittle  ductile  transition. 


2.4  Zone  4:  Tibet 

The  Tibetan  plateau  forms  a  high  topography  region  with  an  elevation  of  4  -  5  km  above  sea 
level  and  is  underlain  by  -70  km  thick  crust.  Earthquakes  in  Zone  4  are  observed  by  both  MKAR 
and  KKAR  at  far-regional  distance,  with  some  overlap  between  the  arrays.  Most  seismicity  is 
located  in  the  Lhasa  Block  (Figure  8),  one  of  the  several  old  suture  zones  that  have  undergone 
north-south  shortening  and  east-west  extrusion  to  form  the  high  topography  and  accommodate 
strain.  It  is  suggested  that  the  entire  plateau  is  underlain  by  cold  lithospheric  mantle  (Priestley  et 
al.,  2006),  however  the  extent  that  the  Indian  plate  subducts  beneath  the  plateau  is  still  debated 
(e.g.,  Curtis  and  Woodhouse,  1997;  Zhou  and  Murphy,  2005).  Low  upper-mantle  S  velocities 
beneath  the  central  and  northern  plateau,  inefficient  Sn  propagation  and  slow  Pn  velocities  are  also 
reported  (e.g.,  Barazangi  and  Ni,  1982;  McNamara  et  al.,  1997). 

Most  seismicity  across  the  plateau  occurs  at  shallow  depths  (<  25  km),  and  deeper  earthquakes 
tend  to  be  poorly  constrained  (Langin  et  al.,  2003).  The  shallow  earthquake  depth  may  be 
indicative  of  a  hot  upper  mantle  below  the  plateau,  creating  a  shallow  brittle/ductile  transition.  A 
seismicity  study  from  the  INDEPTH  III  experiment  (Langin  et  al.,  2003),  shows  that  earthquakes 
predominantly  exhibit  both  normal  and  strike  slip  faulting  mechanisms.  However,  some  do  exhibit 
thrust  mechanisms.  Fault  plane  solutions  are  consistent  with  shortening  in  the  north-south 
direction,  as  well  as  extension  in  the  east-west  direction. 


10 


80‘E 


85'E 


90"E 


95‘E 


100‘E 


105‘E 


40’N 


35N 


30‘N 


75‘ E 


Figure  9.  Zone  4  and  Zone  5.  Red  circles  are  earthquakes  observed  at  far-regional  distance 
by  the  KKAR  array,  blue  circle  are  observed  by  the  MKAR  array,  purple  circles  are 
observed  by  both  arrays.  Gold  circles  are  earthquakes  outside  the  far-regional  range 
from  these  two  arrays.  The  major  suture  zones  across  the  plateau  are  labeled.  The 
Chinese  Lop  Nor  test  site  is  also  labeled  (green  circle). 


2.5  Zone  5:  Eastern  Tarim  Basin 

This  zone  includes  earthquakes  observed  by  the  KKAR  array  at  far-regional  distances  (Figure 
8).  It  includes  part  of  the  Himalayan  Block  of  the  Tibetan  Plateau  and  the  eastern  edge  of  the  Tarim 
basin,  as  well  as  part  of  the  Qaidam  basin.  The  Chinese  Lop  Nor  nuclear  test  site  sits  to  the  north  of 
the  zone.  A  large  left-lateral  fault  system  (Altyn  Tagh)  separates  the  Tibetan  Plateau  from  the 
Tarim  basin  in  this  region  (Hetzel  et  al.,  2004).  Across  the  basins  to  the  plateau  crustal  thickness 
varies  from  48-70  km  (Leech  et  al.,  2010).  Most  earthquakes  have  focal  depths  above  35  km,  with 
a  few  extending  deeper.  Focal  mechanisms  in  this  zone  show  both  left-lateral  strike  slip  motion 
(along  the  boundary  between  the  Tarim  basin  and  the  Tibetan  Plateau)  and  thrust  faulting  on  the 
Tibetan  Plateau. 


2.6  Zone  6:  Western  Mongolia 

Seismicity  in  western  Mongolia  occurs  along  the  Altay  Mountains  (Figure  9).  It  is  recorded  at 
far-regional  distances  by  the  KKAR  array.  The  active  tectonics  of  Mongolia  results  from  the 
collision  of  the  Indian  and  Eurasian  plates.  Most  of  the  active  faulting  and  seismicity  occurs  in  the 
Western  portion  of  Mongolia  along  the  Altay  Mountains.  Within  the  linear  mountain  belt  are 
parallel,  right-lateral  strike-slip  faults  and  some  E-W  thrust  faults.  The  faulting  is  thought  to 
accommodate  shortening  between  western  China  and  Siberia  (Walker  et  al.,  2006). 
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According  to  Bayasgalan  et  al.  (2005),  all  earthquakes  in  Mongolia  have  centroid  depths 
shallower  than  20  km.  In  the  Altay  Mountains,  focal  mechanisms  are  consistent  with  right-lateral 
strike-slip  faulting  on  NW-SE  oriented  faults.  Thrust  mechanisms  occur  on  E-W  faults.  The 
pattern  of  faulting  might  accommodate  vertical  axis  rotation  around  fault  bounded  blocks 
(Baljinnyam  et  al.,  1993).  Seismic  studies  of  the  crustal  structure  of  western  and  central  Mongolia 
have  not  been  yet  been  published.  Results  from  a  2003  seismic  deployment  across  Mongolia  and 
the  Baikal  region  are  currently  in  press.  Most  information  concerning  the  crustal  structure  comes 
from  petrologic  and  thermobarometric  studies  of  young  xenoliths  and  also  gravity  modeling  (e.g., 
Petit  et  al.,  2002).  These  studies  suggest  that  the  crust  is  46  km  under  the  Altay  Mountains. 
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Figure  10.  Map  of  regional  seismicity  for  Mongolia  and  the  Baikal  Rift  Zone.  Earthquakes 
at  far-regional  distance  from  KKAR  (14  - 18°)  are  shown  in  red  and  from  MKAR  in  blue. 
Other  seismicity  is  denoted  with  gold  circles.  Labels  indicate  places  discussed  in  the  text. 


2.7  Zone  7:  Central  Mongolia 

There  is  less  seismicity  in  central  Mongolia  compared  to  the  Altay  Mountains  of  western 
Mongolia.  The  earthquakes  in  central  Mongolia  are  observed  by  the  MKAR  array  at  far-regional 
distances.  They  occur  on  the  eastern  edge  of  the  Hangay  region,  a  Proterozoic  continental  block 
that  is  bounded  on  the  east  by  a  large  right-lateral  strike-slip  fault  system.  The  Hangay  region 
forms  a  topographic  dome  and  is  mostly  aseismic,  although  normal  faulting  events  do  occur  on  its 
northern  edge.  The  crust  beneath  the  Hangay  has  a  maximum  depth  of  50±3  km  (Petit  et  al.,  2002). 
Most  earthquakes  in  the  region  occur  shallower  than  20  km  (Bayasgalan  et  al.,  2005).  Focal 
mechanism  show  right-lateral  strike  slip  faulting  along  north-south  fault  orientations,  consistent 
with  the  local  tectonics. 
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2.8  Zone  8:  Baikal  Rift  Zone 


The  Baikal  rift  forms  an  1 800-km  depression  in  southern  Siberia.  It  has  been  one  of  the  most 
seismically  active  rifts  on  the  planet  for  more  than  the  last  30  Ma,  but  only  exhibits  10-20  km  of 
overall  extension  (0.3  -  0.6  mm/yr).  Geodetic  rates  indicate  much  faster  extension  rates  on  the 
order  of  1  mm/yr  (Suvorov  et  al.,  2002).  Seismicity  here  is  recorded  by  the  MKAR  array  at 
far-regional  distances  (14  -  18°). 

Numerous  studies  have  examined  the  crustal  structure  of  the  rift  and  surrounding  region  (e.g., 
Gao  et  al.,  2004;  Brazier  &  Nyblade,  2003;  ten  Brink  &  Taylor,  2002;  Suvorov  et  al.,  2002; 
Bayasgalan  et  al.,  2005;  Emmerson  et  al.,  2006).  Receiver  function  profiles  reveal  a  crust  that 
ranges  in  thickness  from  37  -  45  km  (Gao  et  al.,  2004).  The  thickest  crust  is  found  to  the  north 
under  the  Siberian  platform  and  south  below  the  Mongolian  fold  and  thrust  belt.  The  thinnest  crust 
is  below  the  central  portion  of  the  Baikal  rift,  where  the  crust  is  35  km  thick.  Upper  mantle  seismic 
velocities  show  a  2  -  5%  slow  anomaly  below  the  rift  axis,  likely  due  to  increased  temperature 
from  extension  related  mantle  upwelling.  The  crust  within  the  rift  and  surrounding  region  has  an 
average  P  velocity  of  6.4  km/s  (Suvorov  et  al.,  2002).  A  detailed  seismic  study  of  the  rift  (ten 
Brink  &  Taylor,  2002)  shows  up  to  4  km  of  slow  sediments  (1-3  km/s)  underlain  by  20  km  of 
with  typical  continental  velocities  (5  -  6.5  km/s)  and  10  km  of  relatively  fast  lower  crust  (7  -  7.4 
km/s). 

Earthquakes  in  this  region  are  mostly  associated  with  extension  and  show  predominantly 
high-angle  normal  faults,  with  the  P  axis  oriented  to  the  NNW.  Some  strike-slip  faulting  is  also 
observed  (Emmerson  et  a\,  2006).  In  the  western  portion  of  the  rift,  earthquakes  occur  above  16 
km  depth,  but  earthquakes  of  up  to  30  km  depth  are  seen  in  the  NE  rift  portion. 


3.  Phase  Characterization  Using  Small-Aperture  Regional  Arrays 

The  regional  seismic  arrays  of  the  International  Monitoring  System  (IMS)  should  be  a  rich  data 
source  for  the  study  of  far-regional  phase  behavior,  since  they  are  all  comprised  of  high-quality 
borehole  seismometers  that  record  high-fidelity  signals.  However,  beyond  regional  distances,  the 
small  apertures  (<  5  km)  and  spatial  configurations  of  these  arrays  limit  their  usefulness  beyond  P 
and  S  first-arrival  onset  picks.  Standard  array  methods  (e.g.,  slant  stacking  and 
frequency-wavenumber  analysis)  cannot  resolve  the  azimuths  and  slownesses  of  primary  and 
secondary  arrivals,  making  confident  arrival  identification  and  classification  difficult. 

The  main  function  of  a  seismic  array  is  to  improve  signal-to-noise  ratio  through  stacking  and  to 
measure  arrival  times  and  directional  parameters  from  the  beam-forming  process.  The 
performance  of  a  particular  array  depends  on  many  factors,  including  site  noise  levels,  signal 
coherence  and  the  array  geometry.  For  optimal  performance  an  array’s  geometry  (aperture  and 
inter-station  spacing)  should  be  tuned  to  the  characteristics  (e.g.,  frequency,  phase  velocity,  and 
direction)  of  the  signals  of  interest.  Array  geometry  tuning  is  typically  a  balancing  act  to 
adequately  address  optimal  cross-array  signal  coherence,  spatial  aliasing,  measurement  resolution, 
as  well  as  the  local  noise  characteristics.  The  trade-offs  between  array  aperture  and  inter-sensor 
spacing  can  be  generalized  in  several  ways.  First,  larger  array  apertures  imply  a  longer 
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triangulation  arm,  which  allows  for  more  confident  back-azimuth  estimates.  A  larger  aperture  also 
decreases  the  likelihood  that  the  local  noise  field  will  be  coherent  across  the  array,  thereby 
improving  signal  gain  levels  achieved  by  stacking.  Arrays  with  larger  apertures  also  produce  more 
precise  measurement  of  signals  with  high  phase  velocities  (steep  incident  angles),  because  the 
travel  time  across  the  array  is  increased.  Conversely,  if  the  array  aperture  is  too  large,  structural 
heterogeneity  across  the  array  and  other  site  effects  can  negatively  impact  signal  coherence 
between  array  elements  and  degrade  array  stacking  performance.  Finally,  large  inter-sensor 
spacing  also  increases  the  possibility  of  wavenumber  aliasing,  since  the  spatial  Nyquist  frequency 
is  related  to  the  maximum  sensor  separation  distance.  Wavenumber  aliasing  increases  the 
uncertainty  associated  with  measurements  of  multiple  arrivals  within  an  analysis  window. 

The  expected  phase  velocities  of  far-regional  P  arrivals  (e.g.,  Pn,  P4 1 0,  P660,  and  depth 
phases)  range  between  8  -  12.5  km/s,  using  the  AK135  reference  model  (Kennett  et  al.,  1995). 
Most  signal  energy  occurs  in  the  1  to  4  Hz  band,  based  on  spectral  analysis  of  our  far-regional 
dataset  at  MKAR  and  KKAR.  This  indicates  that  incoming  signals  have  wavelengths  of  2  -  10  km 
with  corresponding  spatial  frequencies  of  0.5  to  0.1  (units  of  1/km  which  is  equivalent  to  the 
wavenumber/27i).  This  is  illustrated  in  Figure  10,  where  we  show  the  relationship  between 
temporal  frequency  and  spatial  frequency  for  the  expected  range  of  far-regional  phase  velocities. 
Spatial  frequency  is  important  for  understanding  an  array’s  response. 


Figure  11.  Temporal  frequency  versus  spatial  frequency  for  the  expected  range  of 
far-regional  phase  velocities.  The  blue  region  highlights  the  main  frequency  band  of  our 
far-regional  observations.  The  grey  region  marks  the  far-regional  spatial/temporal 
frequency  bands,  which  is  of  interest  to  this  study. 

The  nine-element  borehole  arrays  at  MKAR  and  KKAR  have  similar  designs  (Figures  1 1  and 
12):  both  contain  an  outer  ring  of  five  elements,  an  inner  ring  of  three  elements,  and  one  central 
element.  The  maximum  aperture  at  MKAR  depends  on  the  direction  of  the  incoming  signal.  For 
signals  from  the  southwest/northeast  the  effective  aperture  is  ~3  km  and  for  signals  from  the 
northwest/southeast  the  effective  aperture  is  ~5  km.  The  non-symmetrical  nature  of  the  MKAR 
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array  indicates  that  its  performance  varies  depending  on  the  direction  of  the  incoming  signal.  In 
other  words,  analysis  of  signals  from  the  northwest/southeast  recorded  at  MKAR  will  be 
potentially  better  resolved  in  azimuth  due  to  the  larger  effective  array  aperture,  and  hence  larger 
triangulation  arm.  The  KKAR  array  is  more  symmetric  and  has  an  effective  aperture  of  3  -  4  km, 
depending  on  the  direction  of  the  incoming  signal. 

Examination  of  the  frequency-wavenumber  array  response  functions  for  both  arrays  (Figures 
11  and  12)  shows  that  the  arrays  are  not  optimally  designed  for  full  slowness  and  back- azimuth 
resolution  at  far-regional  distances,  especially  when  considering  multiple  P  arrivals  within  a  5  -  20 
second  analysis  window.  The  main  lobe  of  the  response  functions  has  a  50%  reduction  of  power  at 
the  spatial  frequency  of  0.15  (1/km)  at  MKAR  and  0.2  (1/km)  at  KKAR.  A  simple  rule  of  thumb 
states  that  the  sharper  the  main  lobe,  the  more  precise  the  slowness  measurements.  Thus,  at  the 
temporal  frequencies  of  interest  (1-4  Hz),  the  precision  of  slowness/azimuth  measurements  at 
MKAR  and  KKAR  can  be  negatively  affected,  and  multiple  signal  arrivals  can  become  difficult  to 
discern  with  any  confidence,  especially  when  they  have  similar  slowness  values. 

For  example,  in  the  AK135  reference  model,  an  event  observed  at  15°  epicentral  distance  has  a 
predicted  Pn  arrival  velocity  of  8.15  km/s  and  a  P410  arrival  velocity  of  10.0  km/s.  These 
correspond  to  spatial  frequencies  at  2  Hz  of  0.25  (1/km)  and  0.2  (1/km),  respectively.  Since  the 
main  lobe  of  the  MKAR  response  function  has  a  half-width  of  0. 15,  the  similar  arrival  parameters 
might  smear  the  phases  together,  at  least  for  spectral  array-processing  methods  that  do  not  retain 
arrival-time  information.  At  higher  spatial  frequencies  (shorter  wavelengths)  the  array  response 
functions  exhibit  aliasing  lobes.  Optimal  array  response  functions  have  low-amplitude  side  lobes 
and  a  sharply  peaked  main  lobe.  However,  at  MKAR  an  aliasing  lobe  at  50%  power  occurs  at 
spatial  frequencies  of  -0.45  (1/km)  and  at  KKAR  they  appear  near  0.5  (1/km).  Aliasing  lobes  can 
also  make  it  extremely  difficult  to  accurately  identify  multiple  arrivals,  because  the  high  gain  of 
the  side  lobes  can  amplify  nuisance  signals,  which  can  then  be  misinterpreted.  While 
frequency- wavenumber  response  functions  (as  shown  here)  are  most  directly  applicable  to  spectral 
array  processing  methods,  such  as  the  Capon  method  (Capon,  1969),  in  time-domain  array 
methods  (e.g.,  Nth  root  and  phase-coherence  semblance  stack),  the  effects  of  the  broad  main  peak 
and  side  lobes  manifest  as  smearing  in  the  slowness  and  azimuth  vespagrams,  which  also  reduces 
measurement  precision.  We  will  demonstrate  this  effect  in  the  next  section  of  the  report. 

In  summary,  the  configurations  of  the  MKAR  and  KKAR  arrays  are  not  optimally  designed  for 
analyzing  the  far-regional  P  arrival  suite.  The  primary  difficulty  we  have  encountered  is  the  small 
array  aperture,  which  limits  the  precision  in  the  arrival  measurements  for  the  range  of  far-regional 
phase  velocities  and  frequencies  when  analyzed  with  common  array  methods.  The  lack  of 
precision  makes  confident  phase  identification  difficult,  since  multiple  arrivals  with  similar 
propagation  parameters  occur  within  a  5  -  20  second  arrival  window.  To  mitigate  these 
limitations,  we  present  a  modified  array  processing  scheme  in  the  next  section  (phase-weighted 
semblance  stacking)  that  overcomes  some  of  these  issues.  Conversely,  we  note  that,  in  general,  the 
configuration  of  the  MKAR  and  KKAR  arrays  performs  well  in  regional  waveform  analysis  or  for 
first-arrival  onset  picking  of  teleseismic  phases,  where  straight  stacking  to  reduce  noise  is  the 
primary  objective. 


15 


0 


0.3 

Spatial  frequency  (1/km) 


E 


>s 

O 

c 

o 

=5 

cr 

fi 


CL 

if) 


Easting  (km) 


Spatial  frequency  (1/km) 


-30  -20  -10  0 

Amolitude  resoonse  (dBI 


Figure  12.  MKAR  array  response  function.  Top :  Cross  section  through  response  function, 
showing  the  spatial  frequency  at  which  the  amplitude  gain  is  at  50%  power  loss  and  the 
aliasing  lobe  near  0.45  (1/km).  Bottom  Left :  MKAR  array  configuration;  Bottom  Right : 
MKAR  response  function.  All  amplitudes  are  shown  as  normalized  decibels. 
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Figure  13.  KKAR  array  response  function.  Top :  Cross  section  through  response  function, 
showing  the  spatial  frequency  at  which  the  amplitude  gain  is  at  50%  power  loss  and  the 
aliasing  lobe  near  0.45  (1/km).  Bottom  Left:  KKAR  array  configuration;  Bottom  Right: 
KKAR  response  function.  All  amplitudes  are  shown  as  normalized  decibels. 

4.  Phase-Weighted  Semblance  Stacking  of  Small-Aperture  Array  Observations 

Over  the  course  of  this  effort  we  have  investigated  various  beamforming  methods  and  their 
ability  to  improve  the  analysis  of  far-regional  arrivals  recorded  by  small-aperture  arrays.  The 
investigation  included  both  spectral  frequency-wave  number  methods  (e.g.,  Capon,  1969), 
multiple  signal  characteristic  (MUSIC;  Stoica  and  Nehorai,  1989),  cross-correlation  (Tibuleac  and 
Herrin,  1997),  and  other  slant-stack  beamforming  methods  (e.g.,  Nth  Root,  see  Rost  and  Thomas, 
2002). 


While  frequency-wave  number  methods  (e.g.,  Capon)  can  perform  well  at  measuring  arrival 
slowness  and  back-azimuth,  they  perform  poorly  in  the  case  of  multiple  signal  arrivals  within  a 
single  analysis  window.  The  loss  of  arrival-time  information  on  the  individual  arrivals,  due  to  the 
spectral  formulation,  makes  unraveling  the  P-arrival  suite  nearly  impossible.  At  far-regional 
distances  multi-pathed  arrivals  commonly  interfere  with  each  other  (constructively  and 
destructively)  and  exhibit  similar  arrival  slowness  over  the  5-20  second  duration  of  the  main 
P-wavetrain.  Accurate  arrival-time  information  is  crucial  in  sorting  the  phase  arrival  progression 
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in  such  scenarios.  Windowing  the  arrivals  using  a  Short-Time-Fourier-Transform  method  is  an 
option;  however,  arrival-time  resolution  is  still  problematic  due  to  close  arrival  spacing,  and 
choosing  optimal  analysis  windows  is  difficult  without  affecting  the  quality  of  the  Fourier 
Transform.  In  general,  arrival-time  resolution  is  an  issue  for  most  moving-window  type  analysis 
schemes,  including  slant-stack  Nth-root  and  cross-correlation  implementations.  Methods  such  as 
MUSIC  do  not  suffer  these  limitations  and  can  perform  well  on  far-regional  arrivals;  however,  one 
major  drawback  of  this  method  is  the  need  to  know  a  priori  how  many  arrivals  are  in  the  analysis 
window.  In  many  cases,  emergent  first  arrivals  and  weak  secondary  arrivals  can  be  easily  missed 
without  high-level  scrutiny  by  an  analyst. 

For  the  far-regional  phase  identification  problem,  an  optimal  beamforming  method  needs  to 
enhance  weak  arrivals  and  be  unbiased  in  amplitude  and  phase,  so  that  large  amplitude  arrivals  do 
not  obscure  smaller  amplitude  arrivals.  The  method  should  also  be  able  to  simultaneously 
determine  slowness  and  back  azimuth,  while  preserving  the  phase  arrival  time.  We  have  developed 
a  processing  technique  that  comes  close  to  meeting  these  criteria,  which  we  call  phase-coherence 
semblance  stacking  (PCSS).  This  technique  has  the  signal- to-noise  increasing  power  typical  of 
slant-stack  methods,  but  through  the  incorporation  of  both  a  phase  coherency  measure  and 
amplitude  semblance,  it  produces  a  measure  that  is  highly  sensitive  to  signal  onset  time  and 
overcomes  some  of  the  slowness  and  azimuth  resolution  limitations  of  small-aperture  arrays.  The 
technique’s  sensitivity  to  onset  time  also  holds  for  very  emergent  arrivals. 

In  our  implementation,  both  instantaneous  phase  coherence  and  semblance  are  simultaneously 
determined  for  a  range  of  slowness  and  back-azimuth  combinations  (i.e.,  vector  slownesses), 
resulting  in  a  3-dimensional  (3-D)  phase-coherence  semblance  function  (time,  slowness,  back 
azimuth).  Phase  coherence  and  semblance  are  combined  through  a  linear  mixing,  and  the  two 
coherency  measures  are  complementary.  Instantaneous  phase  coherence  is  an  unbiased  amplitude 
measure  that  is  maximized  for  the  vector  slowness  that  aligns  the  signals  most  closely  in  the 
complex  plane.  Semblance  is  maximized  for  vector  slownesses  at  which  the  total  beam  power 
equals  the  average  beam  power  over  all  traces,  thus  favoring  high-amplitude  correlation  and  high 
signal-to-noise  ratios. 


Instantaneous  phase  is  a  feature  of  complex  signal  analysis,  in  which  the  angular  phase  can  be 
described  for  each  sample  point  in  the  time  domain  through  a  linear  transformation.  To  determine 
the  instantaneous  phase,  we  transform  each  seismic  trace  s(t)  for  station  j  into  its  complex 
representation,  by  combining  the  real  trace  with  its  Hilbert  transform,  such  that 

Xj(t)  =  Sj(t )  +  iH[sj  COl- 

On  the  right  side,  the  term  i  denotes  the  imaginary  unit  and  H[sj  (t)]  is  the  Hilbert  transform  for 
each  seismic  trace  j.  To  equalize  the  power  at  each  sample  point  and  to  limit  Xj(t)  to  the  unit 
circle  in  the  complex  plane,  we  normalize  by  the  magnitude  of  Xj{t), 


Xj(t )  = 


sj(t)  +  iH[sj  (t)] 
sj(t )  +  iH[sj  C0]|" 
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This  limits  the  magnitude  of  Xj(t)  to  values  between  0  and  1.  For  our  purposes,  it  is  more 
convenient  to  rewrite  Xj(t)  using  Euler’s  formula, 

Xj(t)  =  Aj(t)ex p  [i0;(t)] , 

where  each  complex  trace  is  separated  into  its  amplitude  envelope  function,  A(t),  and 
instantaneous  phase  component,  (pit).  Since  we  are  interested  in  optimal  beamforming,  we 
determine  the  phase  coherence  between  all  seismic  traces  by  computing  an  instantaneous  phase 
stack  for  a  range  of  vector  slownesses  @j,  such  that  the  phase  stack,  c,  for  each  time  sample  and  for 
vector  slowness  i,  is  written  as 

M 

Ciid0 1)  =  —  ^  exp  [i6i  (t  +  T{rj,  0i))f- 
7  =  1 

Here  the  term  r(r;-,0j)  is  a  time  delay  that  depends  on  the  array  element  position  vector 
rj(.xk>yk )  and  slowness  vector  dt.  Since  we  are  averaging  the  instantaneous  phase  over  all 
seismic  traces (J  =  1,2,3  ...M),  the  phase  stack,  Cj(0(,  t),  has  a  value  between  0  and  1,  where  0 
means  the  traces  are  completely  out  of  phase  for  that  sample  point  and  slowness  vector,  and  1 
means  the  traces  are  completely  in  phase  and  thus  perfectly  coherent.  The  sharpness  of  the  phase 
stack  can  be  controlled  by  the  exponent  y,  which  can  have  a  value  >  1 .  High  values  can  be  useful 
for  cases  of  low  signal-to-noise  on  some  array  elements,  and  we  have  typically  used  values  of  2  or 
3.  In  some  cases  it  is  also  useful  to  average  the  phase  stack  over  a  time  gate  of  time  samples  m, 
using 

bi+m 

q(0[,  t(n+m/2))  —  —  Ci(Qi,t). 

t=tn 

This  is  particular  useful  to  smooth  the  phase  stack  function  and  in  the  calculation  of  statistics  (e.g., 
the  F-statistic).  In  our  analysis  we  found  that  time  gates  of  5  to  10  samples  produced  good  results. 

A  similar,  but  less  complete  description  of  the  instantaneous  phase  stack  can  be  found  in 
Schimmel  and  Paulssen  (1997).  However,  in  their  implementation,  they  found  it  useful  to  weight 
each  sample  of  the  linear  beam  by  the  phase  stack  at  each  sample  point.  Thus,  in  their 
implementation,  vespagrams  would  be  constructed  by  weighting  the  beam  trace  with  the  phase 
stack.  In  contrast,  we  have  found  it  more  beneficial  to  directly  use  the  phase  stack,  c(t).  This 
permits  a  clearer  picture  of  the  arrival  structure,  since  it  does  not  obscure  small  amplitude  arrivals, 
and  allows  a  better  measurement  of  signal  onset  time.  The  phase  coherence  values  also  lend 
themselves  to  more  straightforward  uncertainty  assessment,  since  a  coherence  value  of  zero 
indicates  that  the  signals  are  completely  out  of  phase  for  a  particular  time,  slowness,  and  back 
azimuth,  while  a  value  of  one  indicates  they  are  completely  in  phase.  The  difference  between  the 
two  vespagram  analysis  methods  is  illustrated  in  Figure  13  for  a  typical  far-regional  event  with 
emergent  first  arrivals  recorded  by  the  MKAR  array  (Table  1). 

Table  1.  ISC  location  for  the  event  shown  in  Figures  13, 14,  and  15. 
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Date 

Origin  Time 

Latitude  (°) 

Longitude  (°) 

Depth  (km) 

nib 

Region 

10/27/2006 

07:55:02 

29.88 

80.04 

10 

4.2 

N.  India 

Figure  13  shows  that  the  vespagrams  constructed  from  weighted  beam  power  are  dominated  by 
the  largest  amplitude  signals  in  the  record.  While  the  emergent  first  arrivals  are  still  apparent 
(between  0-2  seconds  in  Figure  13),  they  appear  significantly  less  energetic  in  the  vespagram 
than  the  dominant  signals,  which  appear  near  the  6-  and  10-second  marks  and  could  easily  be 
missed  in  a  noisier  record.  However,  the  vespagram  slices  representing  the  phase  stack  (i.e.  phase 
coherence)  values  clearly  delineate  the  emergent  first  arrivals  at  the  same  coherence  level  as  the 
larger  amplitude  signals.  It  is  this  aspect  of  the  phase-coherence  method  that  makes  it  particularly 
useful  for  analyzing  far-regional  coda  arrivals. 


Also  apparent  in  Figure  13  is  the  limitation  in  slowness  resolution  of  small-aperture  arrays. 
While  phase-coherence  processing  is  an  improvement  over  other  slant-stacking  methods, 
significant  slowness  smearing  still  occurs,  which  increases  uncertainty  in  phase  identification. 
This  issue  is  slightly  improved  by  including  the  semblance  calculation  (shown  below),  but  is  a 
function  of  the  limited  array  aperture.  It  is  possible  that  further  improvement  will  only  be  realized 
by  increasing  the  aperture  of  the  arrays  and  inter-station  spacing;  for  instance,  by  installing  an 
outer  ring  on  regional  arrays. 


a)  Phase-weight  stack 
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b)  Phase-weight  coherence 
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Figure  14.  Slowness  and  back-azimuth  slices  through  3-D  vespagrams  from  (a)  phase-weighted 
beams  and  (b)  phase-weighted  coherence  values.  The  top  panels  show  the  array  gathers,  the 
middle  panels  show  a  slowness  slice  along  for  a  constant  back-azimuth  of  187°,  and  the  bottom 
panels  show  a  back-azimuth  slice  for  a  slowness  value  of  0.09  s/km.  The  event  was  recorded  by  the 
MKAR  array  at  an  epicentral  distance  of  17°  and  has  a  theoretical  back  azimuth  of  187°  (see 
Table  1). 
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To  further  illustrate  the  differences  between  the  phase-weighted  stack  and  phase  coherence 
vespagrams  we  compare  profiles  (Figure  14)  from  the  back-azimuth  vespagram  slices  shown  in 
Figure  13.  The  profiles  are  taken  along  a  constant  back-azimuth  of  187°  and  represent  beam  power 
or  phase  coherence  as  a  function  of  arrival  time.  In  general,  the  two  profiles  exhibit  a  similar 
shape;  however,  the  profile  of  phase  coherence  shows  much  sharper  onset  times.  This  is 
particularly  true  for  the  emergent  first  arrivals,  as  well  as  later  less  energetic  arrivals.  The 
sharpness  of  the  phase  coherence  function  allows  more  precise  arrival-time  measurements,  since 
the  actual  arrival  onset  is  apparent.  In  the  case  of  beam  power,  arrival  time  is  most  easily  measured 
from  the  peaks  of  the  function;  however,  these  zero-phase  times  would  need  to  be  corrected  before 
being  included  in  a  network-based  event  location  algorithm. 


Phase-weight  stack 
Phase-weight  coherence 


Figure  15.  Profiles  through  the  back-azimuth  vespagram  slices  shown  in  Figure  13.  Profiles  are 
taken  along  a  constant  back-azimuth  of  187°.  The  red  line  represents  the  phase-weighted 
coherence  values,  and  the  blue  line  represents  normalized  phase-weighted  beam  power.  The 
phase-weighted  coherence  profile  shows  much  sharper  signal  onsets. 

While  phase  coherence  stacking  performs  well  in  the  analysis  of  far-regional  arrivals  observed 
on  small-aperture  arrays,  there  is  still  significant  smearing  in  back-azimuth  and  slowness  space. 
This  adds  uncertainty  to  the  array  measurements  and  subsequent  phase  identification.  We  found 
further  improvement  was  possible  by  combining  phase  coherence  and  semblance.  Semblance  is 
defined  by  the  following  equation  for  vector  slowness  9^, 


s”i*dt+rM‘))r 

x,  (t  +  7(r,,e,)) 


The  notation  is  the  same  as  in  the  phase  coherence  stack  description,  with  time  delay  T(rj,  0j)  for 
array  element  position  r;-.  The  semblance  calculation  is  amplitude  based  and  can  be  described  as 
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the  ratio  of  the  total  beam  power  to  the  average  power  across  all  array  elements.  It  most  sensitive  to 
slowness  values  that  produce  high  amplitude  coherence.  It  has  values  that  range  between  0  and  1, 
as  in  the  phase  coherence  stack.  It  is  also  beneficial  to  average  the  semblance  calculation  over  a 
time  gate.  To  combine  the  phase  coherence  stack  and  semblance  we  use  a  linear  mixing  to  produce 
the  phase  coherence  semblance  stack  (PCSS), 


PCSSiiGi.t ) 


Si(di,t )  +  Cj(6>i,  t) 
2 


The  PCSS  is  a  3-D  function  with  parameters  of  slowness,  back  azimuth  and  time.  The 
combined  coherence  values  also  range  between  0  and  1 .  An  example  PCSS  vespagram  is  shown  in 
Figure  15  for  the  same  event  shown  in  Figure  13,  with  similar  slices  through  the  3-D  vespagram. 
However,  the  filter  band  is  slightly  different  (0.5  -  2.85  Hz  in  Figure  13  compared  to  0.8  -  3.5  Hz 
in  Figure  15),  as  are  the  color  scale  and  axis  labels. 


One  of  the  main  advantages  of  combining  the  semblance  with  the  phase  coherence  stack  is  the 
reduction  in  slowness  smearing  and  the  more  clear  separation  of  individual  arrivals.  This  is 
apparent  through  a  comparison  of  Figures  13  and  15.  Both  the  phase  coherence  stack  (Figure  13) 
and  the  PCSS  (Figure  15)  illuminate  the  highly  emergent  first  arrivals,  which  are  hard  to  discern 
with  the  eye.  In  the  PCSS  slowness  vespagram,  two  arrivals  between  0  and  1  seconds  are  apparent 
with  arrival  slownesses  of  -11.5  s/deg  and  -13.0  s/deg.  In  the  straight  phase  coherence  stack 
(Figure  13),  these  appear  as  a  single  arrival.  Individual  arrivals  occurring  in  wave  packets  are  also 
discernible  in  PCSS  slowness  vespagram,  such  as  those  that  arrive  between  1-3  seconds,  which 
exhibit  a  measurable  increase  in  slowness  over  time.  The  arrivals  are  also  apparent  in  the  PCSS 
azimuth  slice,  although  with  a  lower  coherence  value.  This  is  illustrative  of  one  of  the  drawbacks 
of  presenting  the  3-D  vespagrams  as  2-D  slices.  For  example,  the  azimuth  slice  is  determined  for  a 
constant  slowness  value  (or  phase  velocity;  11.1  km/s  in  this  case),  which  is  typically  not  an 
optimum  value  for  each  arrival.  In  our  analyses  we  designed  a  detection  algorithm  to  process  the 
3-D  vespagrams  to  highlight  arrivals  the  exceeded  a  certain  coherence  threshold. 

We  processed  each  event  in  our  far-regional  database  using  the  PCSS  method,  and  the 
resulting  3-D  vespa  grids  were  then  run  through  a  detection  algorithm,  triggering  on  values  that 
exceeded  75%  coherence.  These  detected  arrivals  are  then  stored  in  an  associated  database  to  be 
processed  by  x-p  methods  (discussed  in  a  later  section  of  this  report). 
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Figure  16.  Example  of  phase  coherence  semblance  stack  analysis  using  data  from  the  event 
in  Table  1.  Top:  Slice  through  the  3-D  vespagram  at  a  constant  phase  velocity  of  11.1 
km/s;  Bottom:  Slice  through  the  vespagram  at  a  constant  azimuth  of  187°. 


5.  Wavefield  Continuation  Analysis  to  Determine  Path-Specific  1-D  Velocity 
Models 

The  most  direct  means  of  phase  identification  is  through  matching  arrival  time  and  slowness 
estimates  to  theoretical  predictions.  In  complex  tectonic  regions  such  as  south-central  Asia,  global 
reference  models  perform  poorly  at  far-regional  distances,  and  more  specific  knowledge  of  the 
along-path  and  near-array  earth  structure  is  required  for  confident  primary  and  secondary  phase 
identification.  Therefore,  we  developed  path-specific  1-D  velocity  models  using  two  wavefield 
continuation  analysis  methods  and  the  array  measurements  made  using  the  PCSS  processing  in  the 
previous  section. 
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Wavefield  continuation  processing  decomposes  a  recorded  wave  field  into  its  plane- wave 
components,  where  each  phase  arrival  can  be  described  by  a  function  of  intercept  time  (r)  and  ray 
parameter  ( p )  (e.g.,  McMechan  and  Ottolini,  1980).  In  this  decomposition,  z represents  an  arrival’s 
two-way  travel  time  from  its  turning  point  at  depth  to  a  point  on  the  surface  directly  above  it,  and  p 
is  its  horizontal  slowness.  In  a  spherically  symmetric  earth,  z-p  analysis  can  completely  describe 
ray  behavior. 


A  ray’s  turning  depth  is  directly  related  to  the  earth  structure  and  hypocentral  distance  along 
its  travel  path.  This  is  illustrated  in  Figure  16  for  rays  in  the  iasp91  model  (Kennett  and  Engdahl, 
1991).  While  the  410-km  and  660-km  discontinuities  cause  multipathing  at  epicentral  distances 
less  than  1 8°,  most  rays  from  the  8  -  1 8°  range  sample  above  300  km  depth;  multipathing  produces 
rays  that  sample  multiple  depths.  At  distances  up  to  28°,  most  rays  have  turned  by  700  km  depth, 
thus,  to  sample  the  upper  mantle  transition-zone  structure  requires  events  in  at  least  the  8  -  28° 
range. 


a)  b) 
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Figure  17.  Illustration  of  the  behavior  of  turning  rays  in  the  iasp91  1-D  reference  model,  a) 
Ray  bottoming  depths  expected  for  events  in  the  8°  -  18°  and  18°  -  28°  distance  range 
(red  and  grey,  respectively);  b)  the  corresponding  arrival-time/slowness  (z-p)  curves. 

To  generalize  our  PCSS  array  measurements  and  explore  the  applicability  of  better  1-D  earth 
models,  we  employed  two  separate  z-p  techniques  to  derive  velocity  calibration  profiles.  Both 
techniques  are  based  on  z-p  (i.e.  delay-time/slowness)  methods  that  decompose  wavefields  into 
their  components  parts.  The  first  method  involves  directly  inverting  array-based  delay-time  and 
slowness  measurements  to  solve  for  a  velocity-depth  function  (e.g.,  Carbonell  and  Smithson, 
1994).  The  second  method  applies  wavefield  continuation  techniques  to  array-beam  record 
sections  to  accomplish  the  same  goal  (e.g.,  McMechan,  1984).  While  these  methods  are  not  new, 
they  are  not  typically  applied  to  small-aperture  array  data,  but  rather  to  data  from  long  linear  arrays 
of  seismometers  (e.g  Morozov  et  al.,  2005).  Since  we  are  applying  these  methods  to  records  from 
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multiple  events,  we  must  account  for  the  effects  of  differing  source  parameters  (origin  time,  event 
depth,  and  focal  mechanism),  which  may  be  unknown  or  poorly  constrained.  From  these  exercises 
we  were  able  to  derive  models  for  specific  regions  in  our  study  area  that  are  based  on  our  array 
observations.  We  discuss  our  application  of  the  two  inversion  methods  in  the  next  two  subsections. 


5.1  Direct  Inversion  of  Delay-Time  and  Slowness  Array  Measurements 

The  direct  inversion  method  uses  the  slowness  and  delay-time  measurements  output  by  the 
PCSS  array-processing  as  input  data.  The  delay-time  of  each  arrival  is  converted  to  r  using  the 
epicentral  distance  and  origin  time  information  found  in  the  ISC  catalog.  For  most  events,  the 
PCSS  vespagrams  were  computed  for  the  first  18  seconds  of  the  array  seismograms.  This  window 
length  captures  all  the  P4\  0/P660  arrivals  and  primary  pP  depth  phases,  but  may  miss  pP  phases 
from  the  410-km/660-km  discontinues  and  their  sP  counterparts  for  deep  events.  Figure  17  shows 
an  example  of  the  technique  we  use  to  collect  delay-time  and  slowness  array  measurements  for 
velocity-depth  inversion. 

We  first  groomed  the  outliers,  using  a  set  of  measurements  made  at  MKAR  consisting  of  172 
earthquakes  that  ranged  from  14  -  29°  epicentral  distance.  This  included  earthquakes  extending 
from  the  Hindu  Kush  region  of  Pakistan/ Afghanistan  to  the  Makran  coast  and  Zagros  Mountains 
of  Iran.  The  grooming  removed  everything  outside  the  8.0  -  14.5  (sec/deg)  slowness  range, 
thereby  eliminating  measurements  from  P-to-S  scattered  signals  and  coherent  noise.  We  also 
discarded  measurements  that  exhibited  significant  slowness  smearing  (i.e.,  measurement 
uncertainties)  during  the  PCSS  processing.  Once  the  groomed  measurements  were  collected,  we 
reduce  them  to  a  single  x-p  curve  by  computing  the  error-weighted-mean  r  value  within  evenly 
spaced  slowness  bins  (Figure  17b),  excluding  bins  with  only  a  few  measurements.  Averaging  the  x 
values  within  a  common  slowness  bin  acts  to  correct  for  errors  in  focal  depth  by  accounting  for  the 
±  x  offset  between  P  and  pP.  While  this  is  not  a  perfect  correction  due  to  potentially  missed  phases 
in  the  array  processing,  the  inclusion  of  sP  arrivals,  and  measurement  error,  we  have  found  it  more 
effective  and  feasible  than  correcting  r  on  an  event  basis  using  catalog  depths,  which  typically 
have  large  uncertainties.  The  averaged  x-p  curves  are  the  observations  that  are  input  into  the 
velocity-depth  inversions. 

We  note  that  there  are  some  differences  between  the  averaged  r-p  curve  from  the  array 
measurements  and  the  iasp91  curve.  While  the  two  curves  follow  the  same  general  trend,  the 
averaged  empirical  curve  is  not  nearly  as  smooth  as  the  iasp91  curve.  This  may  reflect 
measurement  error  in  the  array  processing  or  inadequate  averaging,  or  result  from  the  inclusion  of 
events  from  a  region  of  variable  earth  structure.  The  averaged  curve  also  has  inflection  points  that 
are  more  dramatic  and  differ  in  location  from  those  of  the  iasp91  curve.  In  the  averaged  curve  the 
main  inflection  points  occur  near  10  s/deg  and  11  s/deg,  compared  to  9.5  s/deg  (660-km 
discontinuity)  and  1 1.5  s/deg  (410-km  discontinuity)  in  the  iasp91  model.  This  may  be  related  to  a 
shallowing  of  the  mantle  discontinuities  and  a  more  profound  velocity  contrast  across  them.  The 
measured  x-p  data  also  includes  a  break  in  the  curve  near  13.5  s/deg  followed  by  a  reduction  in 
slope.  These  measurements  come  from  events  in  the  14  -  29°  distance  range.  It  is  unclear  what 
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causes  this  break  (e.g.,  measurement  error,  averaging  effect,  or  true  earth  structure),  but  it  will 
significantly  influence  the  shallow  mantle  velocity  gradient. 

a)  b) 


Figure  18.  Sample  r -p  data  from  the  MKAR  array,  a)  The  groomed  data  (blue  circles);  red 
dots  show  the  IASPEI91  (i.e.,  iasp91 )  r-p  curve  for  a  surface-focus  source.  The  range  in  r 
for  a  particular  slowness  is  caused  by  the  difference  in  t  for  earthquakes  at  different 
depths,  rather  than  just  measurement  error,  b)  The  averaged  r-p  curve  from  the  groomed 
data,  evenly  sampled  in  slowness  at  0.15  sec/deg. 

Figure  18  shows  a  preliminary  velocity  profile  determined  from  the  r-p  data  shown  in  Figure 
17.  There  are  many  ways  to  parameterize  this  inverse  problem,  and  the  results  in  Figure  18 
represent  only  one  solution.  To  derive  the  model  we  solved  for  the  thicknesses  of  a  sequence  of 
layers  that  best  fit  the  z-p  observations,  using  a  damped  least-squares  formulation.  For  this 
particular  inversion,  we  linearized  the  problem  by  assigning  each  slowness  p  to  a  layer 
representing  its  maximum  depth  of  penetration  (i.e.  the  ray  turning  point),  which  depends  on  the 
initial  velocity-depth  profile  and  is  updated  after  each  inversion  iteration.  The  under-determined 
system  of  equations  we  invert  is  z  =  G  -h,  where  r  is  a  vector  of  zero-offset  time  observations  and 
h  is  vector  of  model  layer  thicknesses.  Each  row  of  G  is  composed  of  the  vertical  slowness 
through  each  layer  up  to  the  turning  depth  for  each  p  observation.  While  linearizing  the  inversion 
allows  a  direct  estimation  of  model  uncertainties,  the  final  results  tend  to  be  strongly  influenced  by 
the  initial  model. 
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a)  b)  c) 


Velocity  (km/s)  Depth  (km)  Observation  #  (o) 


0.0  0.5  t.O 

resolution 

Figure  19.  r-p  inversion  for  a  velocity  profile  using  the  data  shown  in  Figure  17.  a)  Final 
model  (red)  compared  to  the  starting  model  (black);  b)  the  model  resolution  matrix;  and 
c)  the  data  resolution  matrix. 


5.2  Wavefield  Continuation  Analysis  Using  Array  Beams 

The  method  of  downward  wavefield  continuation  (e.g.,  Clayton  and  McMechan,  1981)  is 
commonly  used  in  exploration  seismology  to  extract  velocities  from  refraction  profile  data.  To 
apply  it  to  passive  source  seismic  data,  the  refraction  profiles  are  replaced  by  earthquake 
seismograms  arranged  according  to  their  distances  and  aligned  on  their  first  arrivals. 

The  wavefield  continuation  consists  of  two  linear  transformations  of  the  data:  1)  slant  stacking, 
which  transforms  the  data  into  a  p-T  space;  and  2)  downward  continuation,  which  involves 
transformation  of  the  p-rau  data  into  a  p-z  plane.  Here  p  is  a  ray  parameter  (horizontal  component 
of  the  slowness  vector),  and  r  is  the  vertical  component  of  the  travel  time  t.  We  perform  slant 
stacking  in  the  time  domain  using  the  following  equation: 

oo 

s(r,  p)  =  j  F(j  +  px,  x)dx , 
o 

where  F  is  a  trace  observed  at  a  distance  x,  and  r  =  t  -  px.  The  p-zau  wavefield  can  be  mapped 
(migrated)  to  the  depth-slowness  (or  depth-velocity)  field  using  the  downward  continuation 
operator  and  a  velocity-depth  function  v(z).  The  downward  continuation  operator  is  given  by: 

'*'(/’>*)  =  }(v(z)2 -p2)l,2dz 

o 

At  regional  and  teleseismic  distances  the  velocities  should  be  modified  using  the  velocity 
flattening  transformation:  v  (z)  =  v(z)  Re  !(Re  -  z)  ,  where  z  is  the  depth  and  Re  is  the  earth  radius. 
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Because  a  velocity  model  has  to  be  chosen  in  order  to  perform  the  downward  continuation,  the 
inversion  procedure  is  done  iteratively.  The  initial  velocity  function(  v,c  )  is  chosen  for  the  first 
iteration,  and  after  the  downward  continuation  is  completed,  a  new  velocity  vf  is  picked  on  the 
maximum  amplitudes  of  the  downward-continued  data.  This  process  is  repeated  using  the  new 
picked  velocity  until  convergence  is  achieved.  For  example  the  convergence  criterion  using  the 
root-mean-square  (RMS)  difference  between  the  initial  and  the  resulting  velocities  becomes 
|vf  -  v.  |  <  s ,  where  e  represents  some  accepted  tolerance  value,  i  is  an  iteration  number,  and  | . . .  | 
denotes  an  Euclidean  norm  operation. 

One  potential  pitfall  of  this  procedure  was  pointed  out  by  Clayton  and  McMechan  (1981).  If 
the  initial  velocity  is  selected  too  high  or  too  low  at  all  depths,  this  process  may  not  converge  to  a 
single  self-consistent  velocity  function.  Instead  it  may  oscillate  between  two  mutually  exclusive 
pseudo-stable  states.  These  two  states  constitute  bounds  on  the  true  velocity  function,  because  one 
of  them  is  too  high,  and  the  other  is  too  low  everywhere.  The  bounds  diverge  as  the  depth 
increases,  since  the  error  accumulates  with  depth. 

To  avoid  oscillatory  behavior,  we  slightly  changed  the  velocity  update  procedure.  Instead  of 
changing  it  to  the  maximum  of  the  envelope,  we  found  the  weighted  average  between  the  velocity 
used  for  the  downward  continuation  and  the  new  velocity  according  to  a  formula: 

V/+i  =«i  v?  +a2vci, 

where  ax  +  a2  =  1,  0  <  a,  <  1,  0  <  a2  <  1  .  We  found  that  this  procedure  dampens  the 
oscillations  and  provides  convergence  to  an  intermediate  stable  state. 

We  applied  the  wave  continuation  technique  to  a  data  “profile”  composed  of  the  earthquake 
observations  at  both  the  KKAR  and  MKAR  arrays.  We  included  shallow  events  spanning  a 
geographic  area  from  Pamir  and  Hindu  Kush  to  the  Afghan  block  and  Eastern  Iran.  This  covers 
epicentral  distances  between  3.4  -  30°  and  a  back-azimuth  range  between  200°  and  250°.  The 
hypocenter  parameters  for  the  data  set  were  taken  from  the  ISC  catalog,  which  was  filtered  to 
exclude  events  with  depths  greater  than  45  km.  Figure  19  shows  the  locations  of  the  events  we 
used  in  this  exercise  (green  and  blue  circles),  and  Figure  20  shows  the  record  section  that  results 
when  the  array  beams  are  sorted  according  to  epicentral  distance  from  the  MKAR  and  KKAR 
arrays.  We  note  that  sparse  seismicity  in  some  segments  of  the  selected  profde  creates  gaps  in  the 
wavefield. 
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Figure  20.  Map  of  events  used  to  construct  a  record  section  (MKAR  events  -  blue  dots 
KKAR  events  -  green  dots).  The  pink  and  red  dots  were  not  used  for  this  exercise. 
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Figure  21.  Record  section  of  the  data  from  the  blue  and  green  events  shown  in  Figure  19. 
The  data  are  plotted  as  a  function  of  the  epicentral  distance  from  the  MKAR  and  KKAR 
arrays. 
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Prior  to  applying  the  wavefield  continuation  technique  we  preprocessed  the  data  using  the 
following  steps: 

a)  Performed  array  beamforming  using  the  back  azimuth  computed  from  ISC  locations  and 
the  velocity  predicted  by  AK135  reference  model  (Kennett  et  al.,  1995); 

b)  Band-pass  filtered  the  array  beams  between  0.4  and  1.0  Hz; 

c)  Applied  static  corrections  to  account  for  poorly  known  event  depths;  and 

d)  Computed  data  envelopes. 

After  we  completed  the  preprocessing,  we  performed  a  slant  stack  of  the  record  section  data,  the 
results  of  which  are  shown  in  Figure  2 1 .  This  slant  stack  represents  the  data  space  for  the  wavefield 
continuation  inversion  process. 
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Figure  22.  Slant  stack  of  the  preprocessed  data  from  the  record  section. 

Figure  22  contrasts  the  results  found  in  the  first  iteration  (left)  versus  those  found  after  the  30th 
iteration  (right).  The  AK135  model  is  used  as  the  starting  input  in  the  first  iteration  (red  line).  The 
crustal  velocity  was  kept  constant  in  the  inversion,  using  values  adapted  from  regional  studies  in 
Kazakhstan  (e.g.,  Kosarev  et  al.,  1993).  In  this  case  the  inversion  was  performed  using  a  strong 
damping  parameter  of  a; =0.1 5.  The  inversion  process  stabilizes  after  about  10-15  iterations,  but 
we  continue  the  iteration  out  to  30  iterations  to  observe  the  behavior  of  the  iterative  process  and 
design  an  appropriate  stopping  criterion  in  the  inversion.  The  results  in  Figure  22  show  that  the 
models  found  at  the  first  and  30th  iteration  do  not  differ  significantly  from  each  other,  mainly  due 
to  the  strong  damping  parameter  that  was  applied. 
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Figure  23.  Velocity-depth  model  results  following  downward  continuation.  In  both  plots  the 
AK135  model  is  shown  as  a  red  line  and  the  dashed  yellow  line  shows  the  final  velocity 
model  picked  along  the  maximum  of  the  transformed  wavefield  amplitude,  a)  Results 
after  a  single  iteration;  b)  results  after  the  30th  iteration.  The  dashed  blue  line  shows  the 
starting  velocity-depth  function  used  for  the  transformation. 

We  show  the  full  set  of  velocity-depth  functions  obtained  for  iterations  10-30  in  Figure  23a. 
Figure  23b  shows  the  root-mean-square  (RMS)  difference  between  the  velocity  used  for  the 
continuation  and  the  reconstructed  velocity.  We  use  the  RMS  difference  as  a  measure  of  change  in 
the  inverse  process,  since  we  do  not  have  a  more  concrete  criterion  to  guide  the  fitting  process.  The 
RMS  difference  between  the  two  velocities  after  the  final  iteration  is  on  the  order  of  0.08  km/s, 
which  indicates  that  a  self-consistent  model  was  achieved.  The  final  velocity-depth  function 
obtained  after  30th  iteration  (red  dashed  line)  shows  clear  410-km  and  660-km  discontinuities.  In 
addition,  there  is  a  noticeable  Lehman  discontinuity,  expressed  as  an  increased  gradient  between 
220  and  240  km. 

To  demonstrate  the  effect  of  the  damping  parameter  we  show  the  results  from  using  ai= 1  (no 
damping)  in  Figures  23c  and  23d.  In  this  case  the  process  oscillates  between  two  pseudo-stable 
solutions.  The  final  RMS  difference  is  in  fact  larger  than  the  RMS  at  the  initial  iteration  (~ 
0.88km/s).  Figure  23c  shows  that  the  resulting  pair  of  the  mutually  exclusively  velocity  functions 
(blue  lines)  indeed  represents  the  bounds  on  the  final  velocity  shown  with  the  red  dashed  line 
(same  as  in  Figure  23a).  However,  the  starting  model  (AK135  model  shown  with  the  grey  dotted 
line)  is  also  bounded  by  the  two  states  at  almost  every  depth.  We  note  that  the  resulting 
velocity-depth  function  has  consistently  lower  velocities  than  AK135  through  the  entire  transition 
zone. 
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Figure  24.  Results  from  applying  wavefield  continuation  to  earthquake  data,  a) 
Velocity-depth  functions  obtained  as  a  result  of  the  iterative  process  for  the  iterations 
10-29  (blue  lines).  The  red  dashed  line  shows  the  velocity  after  the  final  (30th)  iteration, 
and  the  grey  dotted  line  shows  the  AK135  model  for  comparison,  b)  RMS  difference 
between  the  velocity  used  for  the  downward  continuation  and  the  reconstructed  velocity 
at  each  iteration  for  damping  parameter  «/=0.15.  c)  Velocity-depth  functions  found  by 
the  iterative  process,  except  with  damping  parameter  «/=l.  d)  RMS  difference  for 
damping  parameter  «/=l. 

Figure  24  shows  a  comparison  of  our  preliminary  velocity-depth  profile  derived  from  the 
wavefield  continuation  procedure  to  the  IASPEI91  (i.e.,  iasp91 )  model.  Below  400  km  depth,  our 
model  exhibits  moderately  slower  velocities  than  iasp91,  and  a  gradient  zone  at  the  410  km  rather 
than  a  sharp  discontinuity  (Figure  24a).  A  comparison  of  the  slowness  vs.  distance  curves  for  the 
predicted  P  arrivals  for  each  model  are  shown  in  Figure  24b.  Here  the  differences  are  more 
apparent,  and  our  derived  model  shows  more  structure  above  the  410-km  discontinuity  (slowness 
<  ~12  s/deg),  and  a  significant  difference  in  the  caustic  points  of  both  the  410-km  (termination  of 
B-C  branch  in  Figure  24b)  and  660-km  discontinuity  (termination  of  D-E  branch).  Some  of  the 
features  we  observe  may  be  related  to  an  inadequate  distance  sampling  of  our  array  beams,  as  well 
as  the  windowing  procedure  used  in  the  initial  i-p  transformation.  Both  of  these  would  explain  the 
initial  appearance  of  the  predicted  410-km  and  660-km  triplicated  arrivals  at  greater  distances  (18° 
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and  20°,  respectively)  in  our  derived  model  than  the  in  the  iasp91  model,  where  they  appear  at  14° 
and  18°,  respectively. 
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Figure  25.  Comparison  of  the  velocity-depth  profile  from  wavefield  continuation  to  the 
iasp91  reference  model:  (a)  derived  model  (black)  and  IASPEI91  ( iasp91 )  model  (red),  (b) 
Slowness  vs.  distance  comparison  between  the  derived  model  (black)  and  the  iasp91 
model  (red)  for  predicted  P  arrivals. 


6.  Examination  of  Region-Specific  Phase  Behavior  Using  Waveform  Clustering 

To  better  understand  the  far-regional  phase  behavior  not  predicted  by  x-p  ray  methods,  we  used 
the  results  from  our  velocity-depth  inversion  to  synthesize  suites  of  seismograms  that  are  used  as 
part  of  a  waveform  clustering  algorithm.  The  clustering  algorithm  processes  array  beams  to  derive 
'wavefield  templates';  i.e.,  grouped  observations  with  similar  phase  characteristics.  The  wavefield 
templates  are  further  analyzed  in  a  non-linear  fashion  by  comparing  them  with  the  synthetic 
seismograms,  looking  for  quantitative  explanations  for  the  phase  behaviors  not  predicted  by  x-p 
ray  methods. 

Another  motivation  for  this  effort  was  the  strong  similarity  in  waveform  character  that  several 
groups  of  far-regional  earthquakes  in  our  database  exhibited.  This  led  us  to  hypothesize  that  we 
could  group  events  into  semi-discrete  clusters  for  each  region  in  our  study  area  based  on  the 
similarity  of  particular  waveform  characteristics.  By  categorizing  the  dataset  in  this  way,  the  goal 
was  to  develop  a  set  of  template  events  that  generalize  the  phase  behavior  and  phenomenon 
observed  for  each  region. 
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As  an  example  of  the  types  of  waveform  similarities  we  see  in  our  data  set,  Figure  25  (right 
panel)  shows  the  array  beams  (and  cumulative  stack)  for  four  seismograms  from  western  China 
recorded  by  the  MKAR  array.  The  beams  exhibit  similar  waveform  structure  in  both  phase  arrival 
times  and  structure.  This  grouping  of  earthquakes  was  serendipitously  found  in  the  dataset  during 
array  processing,  which  led  us  to  develop  more  rigorous  clustering  technique  to  find  more  such 
groupings. 


Figure  26.  Left  panel  shows  4  earthquakes  recorded  at  MKAR  that  exhibit  similar  phase 
arrival  characteristics.  Right  panel  shows  the  array  beams  (and  cumulative  stack),  which 
exhibit  similar  waveform  structure  in  both  phase  arrival  times  and  structure. 

In  addition  to  the  arrival  times  and  slowness  of  the  P-arrival  suite,  our  phase-coherence 
semblance  stacking  analysis  yields  array  beams,  which  are  suitable  low-noise  records  for  input 
into  waveform  clustering  schemes.  We  explored  ‘partitionaT  clustering  methods  in  which  we 
incorporate  cross  correlation  and  objective  function-based  classification  schemes  to  sort  our 
database  of  array  beams.  In  partitional  methods,  each  waveform  is  evaluated  for  membership  to 
particular  cluster  groups  by  its  degree  of  similarity  to  other  cluster  members  in  the  group.  The 
clustering  method  we  favored  is  referred  to  as  fuzzy-clustering  (Bezdek,  1981).  In  our 
implementation,  fuzzy-clustering  finds  similar  seismograms  based  on  arrival  alignment  by 
minimizing  the  objective  function 

Jm  =  ^  ^  Uy  |xi  —  42 

i= 1  ./=! 


In  the  above  function,  i/”'  defines  the  membership  criteria  for  each  cluster  group,  x,  are  the 

waveforms  to  be  clustered,  and  Cj  are  the  cluster  centers  (or  template  array-beam  events) 
determined  from  the  set  of  cluster  members.  The  cluster  centers  are  determined  by  computing, 

Cj  =  ^ - , 

TK 

i= i 

and  the  membership  values  are  determined  from 


34 


1 


uu  = 


c 

Z 

k= 1 


X,.  —  C  ■ 


\/(m- 1) 


■ ,  m  >  1, 


where  m  is  the  fuzziness  value.  Larger  values  of  m  relax  the  criteria  for  membership  to  a  particular 
cluster  group.  The  cluster  centers  change  as  new  members  are  added  and  subtracted  to  the  group, 
and  limited  control  over  the  clustering  results  can  be  made  by  predefining  the  cluster  centers.  For 
each  waveform  and  each  cluster,  the  membership  function  has  a  value  between  0  and  1 ,  where 
values  closer  to  1  reflect  a  stronger  similarity  of  the  waveform  to  the  cluster  center. 

The  clustering  is  carried  out  as  an  iterative  optimization  of  the  objective  function  Jm,  in  which 
the  cluster  centers  (cj)  and  membership  matrix  ( u'J  )  are  updated  at  each  iteration  step.  Iterations 

are  stopped  for  max(/ |m^+i  -  w*|j<  s,  where  s  is  the  analyst-defined  convergence  limit  (between 
0  and  1)  and  k  is  the  iteration  step. 

Our  general  algorithm  for  determining  waveform  clusters  using  the  above  scheme  is  composed 
of  the  following  steps: 

1)  Choose  a  set  of  waveforms  to  be  clustered  and  the  number  of  clusters  to  search  for. 

2)  Cross  correlate  the  waveforms  so  that  they  are  time  aligned. 

3)  Initialize  the  membership  matrix  (u™,k  =  0 ). 

4)  Calculate  the  cluster  centers  (  ck  )  for  the  current  membership  matrix. 

5)  Update  the  membership  matrix  (  n"‘  k ,k  -  k  + 1 )  with  the  current  cluster  centers. 

6)  Calculate  the  convergence  error. 

7)  Stop  the  iterations  if  the  convergence  limit  is  reached,  otherwise  return  to  step  4. 

Figure  26  shows  an  example  cluster  group  determined  using  a  set  of  167  array  beams  from 
varying  geographic  regions.  These  array  beams  were  clustered  into  25  separate  groups,  and  each 
group  was  then  further  separated  into  2-5  groups.  In  this  particular  example,  15  seconds  of  each 
waveform  were  processed,  beginning  0.5  seconds  prior  to  the  initial  arrival  and  including  14.5 
seconds  of  P  coda.  Before  clustering,  the  waveforms  were  time-aligned  to  the  maximum 
correlation  lag  in  this  15  second  window.  For  many  of  the  seismograms,  this  resulted  in  the 
waveforms  being  aligned  along  the  P410  arrivals,  since  these  arrivals  typically  have  the  largest 
amplitude  in  the  window.  In  latter  analyses,  we  limited  the  cross-correlation  to  the  first  arrival 
window  to  better  understand  the  effects  of  the  410  discontinuity  on  arrival-times  and  depth  phases. 

As  evidenced  by  the  members  of  the  sub-group  in  Figure  26,  the  clustering  algorithm  grouped 
waveforms  that  generally  have  similar  P410  phase  structure,  regardless  of  the  P410  arrival  time 
(arrival  time  relative  to  the  first  arrival),  and  that  tend  to  have  emergent  first  arrivals.  When  sorted 
by  increasing  correlation  lag  time,  the  group  members  show  a  positive  move-out  time  of  the  P410 
phase.  For  a  common  earthquake  depth,  a  positive  move-out  of  the  P410  arrival  is  expected  for 
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decreasing  earthquake  distance,  but  for  this  group,  earthquake  distance  does  not  correlate  with  the 
move-out  time  (Figure  26:  right  column).  There  are  several  factors  that  could  explain  this  oddity. 
Since  the  waveforms  grouped  in  Figure  26  are  from  different  regions,  variations  in  the  depth  of  the 
410-km  discontinuity  would  affect  the  arrival  time  of  the  P4 10  phase.  Increasing  earthquake  depth 
would  also  affect  the  P410  arrival  time.  A  more  mundane  explanation  may  be  that  in  some  instance 
the  large  amplitude  phase  in  Figure  26  is  not  the  P410  phase. 

P410 


time 


Figure  27.  An  example  cluster  group  found  from  processing  167  array  beams.  The 
waveforms  exhibit  similar  P410  phase  structure  and  are  time-shifted  based  on  the 
maximum  correlation  lag  time.  Red  dots  mark  the  first  arrival,  and  numbers  on  the  left 
column  are  the  event  ids.  On  the  right  is  shown  the  earthquake  distance  (Delta;  degrees), 
the  station-to-event  azimuth  (SEAZ;  degrees),  and  the  cross  correlation  lag  time  (Tlag; 
seconds). 
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While  the  cluster  group  shown  in  Figure  26  is  interesting,  we  wanted  to  find  cluster  groups  that 
exhibit  both  a  similar  arrival-time  and  phase  structure.  We  noted  that  the  cross-correlation  of  long 
time  windows  tended  to  align  the  waveforms  along  P410  phase.  This  biased  the  clustering  results 
towards  groups  with  similar  P410  phase  structure,  to  the  detriment  of  other  phases  in  the  records. 
In  order  to  minimize  this,  we  first  modified  our  algorithm  to  time  align  on  the  first  arrival,  rather 
than  the  largest  amplitude  arrival.  Second,  we  explored  different  minimization  measures.  The 
cluster  group  in  Figure  26  was  found  using  an  L2  norm,  but  this  measure  can  be  sensitive  to  noise 
in  the  data  and  variances  in  arrival  amplitude  (and  total  power)  between  cluster  members.  The 
variance  in  arrival  amplitude  can  be  a  result  of  different  focal  mechanisms  for  the  different  cluster 
members.  We  found  improved  results  by  optimizing  the  objective  function  with  LI  norm,  which  is 
less  sensitive  to  data  noise,  and  a  semblance  measure,  which  is  less  sensitive  to  amplitude  variance 
and  polarity  swings.  Our  final  implementation  tracked  four  measures  to  compare  the  trade-offs  of 
each:  LI  norm,  L2  norm,  semblance,  and  the  maximum  cross-correlation  coefficient. 

Other  modifications  we  have  made  to  our  algorithm  include  processing  the  Hilbert  envelopes 
of  the  waveforms  (calculated  after  cross-correlation),  and  processing  waveforms  from  specific 
regions  rather  than  the  whole  dataset  at  once.  Using  the  Hilbert  envelopes  removes  the  amplitude 
polarity  information  from  the  waveforms,  which  allows  earthquakes  with  difference  source 
mechanisms  but  similar  arrival  time  structure  to  be  grouped.  This  is  an  important  factor  when 
attempting  to  group  events  from  specific  geographic  regions  where  the  source  mechanisms  where 
highly  varied  for  small  magnitude  events. 

For  example,  Figure  27  shows  the  cluster  group  and  associated  wavefield  template  that  results 
from  applying  our  modified  waveform  clustering  algorithm  to  a  set  of  events  recorded  at  the 
MKAR  array  in  the  14—18°  distance  range.  In  this  example  the  algorithm  clusters  78  events  from 
northern  Pakistan  into  5  groups,  with  12-22  members  per  group.  To  account  for  some  source 
effects,  the  cluster  membership  is  applied  to  the  Hilbert  envelope  of  the  beams,  although  we 
plotted  the  raw  beams  in  Figure  27.  We  chose  the  cluster  groups  based  on  an  Ll-norm  distance 
measure  and  aligned  the  array  beams  prior  to  clustering  using  multi-channel  cross  correlation 
(VanDecar  and  Crosson,  1990)  applied  to  the  first  four  seconds  of  each  beam  signal.  As  previously 
mentioned,  we  found  the  alignment  of  the  beams  to  be  beneficial  to  the  clustering  procedure. 
However,  the  length  of  the  cross-correlation  window  is  an  important  variable.  For  example,  if  the 
window  is  too  long,  the  array  beams  align  on  the  maximum  amplitude  signal  in  the  record.  For 
some  beams  the  maximum  amplitude  arrival  may  occur  early  in  the  seismogram,  while  for  others  it 
occurs  much  later.  This  can  result  in  some  unusual  clustering  of  events  that  seems 
counter-intuitive.  For  this  reason,  we  find  clustering  methods  based  strictly  on  cross  correlation,  as 
used  in  other  types  of  studies  (e.g.,  Menke  1999;  Ferretti  et  al.,  2005;  Barani  et  al.,  2007)  were  not 
suitable  for  a  dataset  such  as  ours. 
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Figure  28.  Wavefield  clustering  to  generate  a  wavefield  template.  The  bottom  panel  shows 
the  cluster  members  sorted  by  epicentral  distance  (shown  on  the  left).  The  green  band 
highlights  the  region  of  cross  correlation.  The  top  panel  shows  the  computed  wavefield 
template  (black)  and  1  -a  deviation  (red). 

In  the  example  shown  in  Figure  27,  the  resulting  cluster  members  include  earthquakes  that  are 
near  the  14.3°  distance  range,  but  some  more  distant  earthquakes  are  also  included.  The 
earthquake  depth  range  is  also  variable,  but  may  be  due  to  poorly  constrained  catalog  depths.  In 
general,  the  cluster  members  show  consistent  waveform  structure  in  the  cross-correlated  portion  of 
the  signal.  Later  arrivals  show  more  variability  between  the  cluster  members,  but  still  appear  to 
align  well  at  approximately  8  seconds  and  later.  The  most  variability  occurs  between  4  and  8 
seconds  delay  time,  where  some  records  show  large  amplitude  arrivals  while  other  show  little 
signal.  Some  of  this  variability  is  likely  due  to  differences  in  earthquake  distance,  triplicated 
arrivals,  and  the  presence  of  depth  phases. 

A  further  example,  shown  in  Figure  28,  depicts  some  of  the  template  beams  that  result  from 
other  cluster  analysis  runs.  This  particular  example  shows  the  results  from  clustering  80 
earthquakes  within  a  back- azimuth  range  of  190  -  230°  from  the  MKAR  array,  which  includes 
events  from  northern  India  and  northeast  Pakistan.  The  clustering  produces  4  template  beams, 
which  qualitatively  appear  to  exhibit  similar  polarity  and  arrival  structure  for  the  first  4  seconds 
signal.  This  likely  implies  a  similar  source  mechanism  for  most  of  the  80  earthquakes. 
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Figure  29.  Template  beams  from  the  cluster  analysis  of  80  MKAR  recorded  events  in  the  190 
-  230°  back-azimuth  range.  The  top  4  traces  are  the  individual  template  beams,  which  are 
shown  overlain  at  the  bottom  of  the  panel.  The  traces  are  aligned  on  the  Pn  phase  and 
filtered  in  the  0.03  -  1.0  Hz  band.  The  P410  phase  is  also  labeled.  At  Trace  4  (red)  P410 
arrives  ~10  seconds  after  Pn,  while  at  Trace  1  (black)  P410  arrives  ~6  seconds  after  Pn. 

After  the  first  4  seconds,  the  template  beams  are  more  dissimilar,  which  is  likely  related  to 
differences  in  earthquake  depth  and  distance  between  the  events  that  compose  each  template 
beam.  The  main  feature  in  this  part  of  the  record  is  the  large-amplitude  P4 1 0  wave  packet  that 
arrives  between  6  and  1 1  seconds  after  Pn.  The  arrival-time  moveout  of  this  wave  packet  is  related 
to  earthquake  distance,  and  this  feature  should  be  verifiable  by  modeling.  In  general,  the  clustering 
analysis  appears  to  produces  reasonable  clustering  groups.  However,  our  main  objective  was  the 
resulting  template  beams,  which  we  view  as  generalized  waveforms  that  are  representative  of 
particular  source  regions. 

A  primary  goal  is  to  be  able  to  explain  of  the  waveform  structure  of  the  template  beams  as  a 
function  of  earthquake  depth,  earthquake  distance,  and  the  depth  of  the  410-km  discontinuity.  To 
fit  the  template  beams  we  devised  a  grid-search  scheme  that  finds  the  minimum  L2- norm  between 
the  individual  template  beams  and  a  suite  of  synthesized  waveforms.  We  used  a  reflectivity 
method  to  produce  the  synthetic  waveforms  (Kennett,  1988).  The  synthetics  were  generated  for  a 
distance  range  of  14  -  1 8.8°  at  an  interval  of  0. 1 0  and  an  earthquake  depth  range  of  2  to  60  km  at  an 
interval  of  2  km.  For  a  single  source  mechanism,  this  modeling  exercise  produces  1470  synthetic 
waveforms.  Figure  29  shows  an  example  earth-flattened  velocity  model  with  a  70-km  thick  crust, 
which  tapers  into  the  AK135  reference  model  in  the  upper  mantle.  The  model  also  includes  a 
gradient  zone  on  the  410/660  km  discontinuities,  thereby  incorporating  the  results  of  the  r-p 
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transformations  developed  in  the  previous  section.  Some  example  synthetic  waveforms  from  this 
model  are  shown  in  Figure  30. 
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Figure  30.  Velocity  model  used  to  generate  synthetic  waveforms.  This  is  an  earth-flattened 
model  that  tapers  into  AK135  in  the  upper  mantle.  P  velocity  is  denoted  with  a  solid  line, 
and  S  velocity  with  a  dashed  line. 

The  grid-search  minimization  algorithm  filters  the  synthetic  waveforms  to  the  same  frequency 
band  as  the  template  beam.  The  synthetics  are  then  cross-correlated  to  the  template  beam  to 
time-align  the  waveforms,  using  a  1  -  4  second  window  after  the  first  arrival.  An  L2- norm  is 
computed  for  each  synthetic  using  a  12  - 18  second  time  window  and  the  best- fitting  solutions  are 
found. 

Results  from  the  grid-search  minimization  algorithm  are  shown  in  Figures  31  and  32,  for  two 
template  beams  shown  in  Figure  28  (beam  templates  2  and  3).  In  our  current  implementation  the 
only  free  parameters  are  earthquake  depth  and  distance,  and  our  velocity  model  is  the  one  shown  in 
Figure  29.  The  source  mechanism  is  a  north-south  striking  thrust  fault  with  a  45°  dip.  While  this  is 
a  common  mechanism  observed  for  northeast  Pakistan,  this  particular  region  also  exhibits 
strike-slip  faulting,  normal  faulting  and  oblique-slip  faulting. 
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Figure  31.  Example  synthetic  seismograms  for  the  waveform  template-fitting  grid  search. 
This  example  shows  a  seismogram  gather  for  a  range  of  event  depths  but  constant 
epicentral  distance.  Note  the  moveout  of  the  P410  phases  as  a  function  of  earthquake 
depth,  and  the  interference  between  depth  phases  and  P410  for  events  depths  between  10 
and  20  km. 

Figure  31  shows  the  grid-search  results  for  template  beam  number  2  (see  Figure  28).  In 
general,  the  best- fitting  synthetic  waveform  captures  the  basic  structure  of  the  template  beam  (top 
panel  of  Figure  31).  The  first  arrival  and  the  .P410  arrival  packets  fit  well,  although  the  synthetic 
does  not  precisely  match  all  the  features  of  the  first  arrival.  Smaller  amplitude  arrivals  are  more 
poorly  fit,  but  do  exhibit  similar  structure  as  the  template  beam.  In  the  bottom  panel  of  Figure  31, 
multiple  minima  are  observed,  indicating  that  multiple  synthetics  fit  the  data  equally  well.  There  is 
a  range  of  distance  solutions  for  an  earthquake  depth  of  22  km  and  another  for  an  earthquake  depth 
of  28  km.  Multiple  minima  are  not  surprising  given  the  nature  of  the  grid  search. 
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Figure  32.  Example  of  the  waveform  fitting  for  template  beam  number  2  shown  in  Figure  28. 
The  top  panel  shows  the  template  beam  (black)  and  the  best-fit  synthetic  (red  dashed). 
The  bottom  panel  shows  the  root-mean-square  residual  for  each  synthetic  waveform  as  a 
function  of  variable  earthquake  depth  and  distance.  The  warm  colors  represent 
waveforms  that  fit  the  data  better. 

Figure  32  shows  the  grid  search  results  for  template  beam  number  3  (Figure  28).  In  this  case, 
the  best-fitting  synthetic  also  captures  the  general  structure  of  the  template  beam.  However,  the 
overall  fit  is  worse  than  the  example  of  Figure  3 1 .  While  the  timing  of  first  arrival  and  P4 1 0  arrival 
is  well  matched,  their  amplitude  is  not.  Some  of  the  amplitude  mismatch  may  be  caused  by  an 
inappropriate  choice  for  the  source  mechanisms.  The  consequence  of  the  poorer  waveform  fit  is 
multiple  solutions,  which,  for  this  example,  include  earthquake  depths  that  range  from  26  to  60  km 
and  a  distance  range  of  14.2  to  18°. 


42 


1 — 1 — 1 — ^ — i — 1 — 1 — 1 — I — 1 — 1 — 1 — r- 

0  4  8  12 

time  (s) 


4  L 


SO  4 — i - I - I - 1 - * - ’ - I - r- 

14  16  18 

distance  (deg) 


l — ' - 1 - ' - 1 — 

1.0  1.2  1.4  1.6 

RMS  residual 

Figure  33.  Example  of  waveform  fitting  for  template  beam  number  3  shown  in  Figure  28. 
The  top  panel  shows  the  template  beam  (black)  and  the  best-fit  synthetic  (red  dashed). 
Bottom  panel  shows  the  root-mean-square  residual  for  each  synthetic  waveform  where 
earthquake  depth  and  distance  vary.  The  warm  colors  represent  waveforms  that  fit  the 
data  better. 

In  general  the  quality  of  the  waveform  fits  were  quite  surprising  given  the  crudeness  of  the  grid 
search  algorithm.  This  suggested  that  the  waveform  templating  scheme  was  able  to  capture  the 
generalized  far-regional  arrival  structure  for  particular  regions.  One  element  of  the  grid  search  that 
remains  to  be  addressed  is  the  best  way  to  select  a  source  mechanism  for  the  synthetics.  The 
waveform  templates  are  essentially  weighted  stacks  of  groups  of  seismograms  that  likely  have 
differing  source  mechanisms,  and  it  is  unclear  how  these  differing  mechanisms  are  reflected  in  the 
templates.  A  more  complicated  approach  would  be  to  try  a  source  inversion  prior  to  performing  the 
grid  search. 
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7.  Studies  of  Near- Array  Earth  Structure 

In  a  supplemental  task  undertaken  during  this  project,  we  studied  the  effects  of  near-receiver 
structure  on  the  wavefield  to  separate  it  from  path-related  effects  that  arise  from  interaction  with 
the  mantle  transition  zone.  To  accomplish  this  we  applied  teleseismic  array  calibration  techniques, 
including  receiver  functions,  polarization  analyses  and  analysis  of  PcP  arrivals,  using  earthquake 
data  observed  30  -  90°  from  the  MKAR  and  KKAR  arrays.  The  purpose  of  using  teleseismic 
signals  for  these  analyses  was  to  limit  the  effects  from  the  interaction  with  the  mantle  transition 
zone  and  produce  array  calibration  information  that  is  independent  of  the  far-regional 
observations. 

We  collected  3-component  broadband  data  from  over  300  teleseismic  events  for  each  of  the 
arrays  and  performed  polarization  analysis  to  determine  a  back-azimuth  residual.  The  residuals 
represent  the  deviations  from  the  great  circle  path  of  the  published  epicenter.  The  back-azimuth 
residuals  are  plotted  as  a  function  of  event  position  in  Figure  33  for  both  MKAR  and  KKAR.  The 
residuals  have  a  range  of  ±20°,  and  a  fair  amount  of  scatter.  As  Figure  33  demonstrates,  there  is  a 
paucity  of  information  in  at  least  one  of  the  quadrants  surrounding  the  arrays,  but  there  are  enough 
data  to  postulate  a  dipping  Moho  at  KKAR  that  causes  a  predominantly  negative  residual  to  the 
southeast.  This  would  be  sensible  in  terms  of  the  geology  near  the  array,  which  is  near  the  edge  of 
deepening  crust.  The  results  at  MKAR  are  less  clear-cut,  as  the  residuals  vary  significantly  as  a 
function  of  azimuth.  This  may  indicate  a  variety  of  conditions  beneath  the  array,  including  a  flat 
Moho  with  significant  local  crustal  heterogeneity.  We  have  found  that  the  far-regional  data  also 
indicate  significant  back-azimuth  residuals  at  MKAR.  One  of  the  objectives  of  this  analysis  was  to 
determine  a  correction  factor  for  the  near-array  structure  (Moho  dip)  that  could  be  used  to  improve 
the  far-regional  array  processing  measurements  of  slowness  and  back  azimuth.  However,  common 
methods  (e.g.,  Niazi,  1966)  applied  to  determine  the  dip  of  the  Moho  below  the  arrays  with  these 
data  resulted  in  highly  non-unique  solutions,  due  to  paucity  of  data  and  outliers.  Assuming  a  single 
dipping  interface,  for  KKAR  we  get  a  strike  solution  of  N45°E  and  a  dip  of  9  ±  6°.  For  MKAR  we 
were  unable  to  determine  a  solution,  which  led  us  to  examine  other  methods  to  determine  the 
near-array  structure,  which  we  describe  in  the  following  subsections. 
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Figure  34:  Back-azimuth  residuals  found  from  polarization  analysis  using  a  large  data  set  of 
teleseismic  events  observed  at  MKAR  (right)  and  KKAR  (left).  The  dashed  line  marks  the 
far-regional  distance  limit.  The  residual  pattern  at  KKAR  is  more  coherent,  which  likely 
indicates  dipping  Moho  structure  beneath  the  array. 


7.1  Near-Array  Earth  Structure  from  Receiver  Functions 

Earth  structure  below  an  array  can  have  a  profound  effect  on  the  observed  back-azimuth  and 
slowness  measurements.  Near-array  structure  is  typically  accounted  for  by  correcting  the 
array-based  measurements  to  theoretical  values.  For  example,  published  earthquake  epicenters  are 
used  to  determine  a  theoretical  back- azimuth  angle  to  the  array,  and  a  1-D  Earth  model  can  be  used 
to  compute  expected  ray  parameters.  In  the  case  of  a  single  dipping  interface  (e.g.,  the  Moho),  the 
azimuth  and  slowness  residuals  relative  to  the  theoretical  values  have  a  sinusoidal  pattern  with 
respect  to  back  azimuth.  The  amplitude  and  phase  of  these  patterns  can  be  used  to  solve  for  dip  and 
strike  of  the  interface,  and  a  correction  can  be  computed  for  each  event  based  on  its  observed 
back-azimuth  angle.  This  method  works  well  when  the  orientation  of  the  interface  can  be  uniquely 
determined,  if  two  conditions  are  met:  1)  the  near-array  structure  must  be  the  phenomenon  causing 
the  anomalous  residuals  (i.e.,  there  are  no  other  path-specific  effects);  2)  the  residuals  must  exhibit 
systematic  behavior  with  respect  to  back  azimuth. 

We  observe  both  back- azimuth  and  slowness  anomalies  at  the  KKAR  and  MKAR  arrays, 
where  azimuth  residuals  are  20°  in  some  instances  (Figure  33).  However,  most  far-regional  events 
occur  to  the  south  of  the  arrays  with  no  events  to  the  north.  This  makes  typical  array  calibration 
difficult  and  highly  uncertain.  It  is  also  not  clear  how  many  of  the  observed  anomalies  are  due  to 
near-array  structure.  To  try  to  obtain  a  more  accurate  accounting  of  the  near-array  structure,  and 
thus  improve  our  array  measurements  with  a  correction,  we  employed  receiver-function 
techniques  to  image  the  structure  below  the  arrays.  Our  goal  was  to  separate  near-array  effects  on 
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the  azimuth  and  slowness  measurements  from  the  along-path  effects  by  accounting  for  the 
near-array  structure  independent  of  the  array  measurements. Determining  near-array  structure  from 
receiver  functions  has  potentially  several  benefits  over  array-based  methods:  1)  only  the  near-array 
structure  is  imaged,  and  far  off  path  effects  have  no  influence,  2)  earth  structure  is  not  determined 
from  the  measurements  to  be  corrected  (the  array  measurements  in  our  case),  and  3)  teleseismic 
events  used  in  constructing  the  receiver  functions  are  more  azimuthally  distributed  compared  to 
the  far-regional  events. 


Figure  34  shows  the  effects  that  dipping  structure  has  on  receiver  function  record  sections.  In 
this  synthetic  case,  the  earth  model  is  a  simple  dipping  Moho  (strike  =  0°  and  dip  =  10°)  over  a  half 
space.  The  receiver  functions  are  computed  for  an  azimuth  range  of  0  -  360°  and  a  constant  ray 
parameter. 
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Figure  35.  Synthetic  receiver  functions  images  for  a  10°  dipping  Moho  over  a  half-space. 
Left:  tangential  receiver  functions;  Right:  Radial  receiver  functions.  Positive  amplitudes 
are  red  and  negative  amplitudes  are  blue.  The  receiver  functions  are  aligned  on  the 
zero-phase  onset  time,  and  the  main  phases  are  labeled. 

As  seen  in  Figure  34,  the  dipping  Moho  structure  causes  both  amplitude  and  arrival-time 
effects.  The  main  Ps  phase  from  the  Moho  shows  an  amplitude  decay  in  the  radial  components  (SV 
energy)  and  a  varying  amplitude  polarity  (sinusoidal)  in  the  tangential  components  (SH  energy). 
The  reverberations  mimic  the  polarity  reversals,  but  also  exhibit  a  time  delay  that  depends  on 
whether  the  signal  is  traveling  up-dip  or  down-dip.  The  phase  and  amplitude  of  the  sinusoidal 
polarity  changes  are  a  function  of  the  strike  and  dip  of  the  interface  and  thus  deterministic. 

For  the  MKAR  and  KKAR  arrays,  we  computed  receiver  functions  using  a  time-domain 
deconvolution  with  data  from  the  single  3-component  instrument  installed  at  each  array.  Figure  35 
shows  distribution  of  teleseismic  events,  displayed  as  the  Moho  Ps  conversion  points,  around  each 
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array.  Several  hundred  receiver  functions  are  computed  for  each  array.  The  distribution  for  the 
MKAR  array  is  more  azimuthally  complete  than  the  KKAR  array. 


MKAR 


KKAR 


Figure  36.  Azimuthal  distribution  of  teleseismic  events  used  in  the  receiver  function 
calculations.  The  green  dots  mark  the  Moho  conversion  point  for  each  of  the  receiver 
functions. 

We  compute  radial  and  tangential  receiver  functions  using  a  time-domain  deconvolution 
method,  and  then  form  azimuthal  record  sections  from  the  receiver  functions.  The  record  sections 
(Figure  36  and  37)  are  constructed  by  stacking  receiver  functions  at  constant  azimuth  spacing.  To 
account  for  variation  in  ray  parameter  of  the  teleseismic  events  and  to  minimize  scattered  energy, 
we  use  a  2-D  Gaussian  weighting  function,  where  the  weight  is  a  function  of  proximity  to  the 
back-azimuth  stacking  point  and  ray  parameter.  In  this  manner,  we  are  able  to  construct  evenly 
sampled  (in  back  azimuth),  smoothed  record  sections  that  can  be  used  in  subsequent  analysis. 

The  main  Ps  converted  signal  at  both  MKAR  and  KKAR  has  a  delay  of  ~3.5  seconds  for  the 
high-frequency  band.  It  is  most  obvious  in  the  radial  component  sections.  At  different  back 
azimuths  it  shows  some  variation  in  arrival  time  and  amplitude,  indicating  a  non-planar  conversion 
interface.  Preceding  this  phase  is  a  large  amplitude  negative  pulse,  which  may  be  related  to  a 
low-velocity  zone  near  the  surface.  Later  arriving  reverberations  are  also  apparent  in  both  the 
MKAR  and  KKAR  sections  and  show  variations  with  back-azimuth.  Coherent  energy  is  apparent 
on  the  tangential  component  sections,  which  is  typically  indicative  of  dipping  structure  or,  in  some 
cases,  anisotropy.  Coherent  energy  on  the  tangential  components  is  clearest  in  KKAR  record 
sections,  where  amplitude  polarity  changes  are  also  seen.  However,  observed  phases  and  the 
coherency  of  the  record  sections  seem  to  have  frequency  dependence.  At  lower  frequencies,  the 
main  Ps  phase  is  less  prominent  and  the  tangential  record  sections  are  less  coherent  (not  shown). 
This  may  be  related  to  the  nature  of  the  velocity  contrasts  across  the  discontinuities  below  the 
arrays,  and  will  likely  not  be  sorted  out  from  our  straightforward  analysis. 
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Figure  37.  MKAR  receiver-function  record  sections.  The  receiver  functions  are  aligned  on 
the  minimum-phase  time  and  have  been  filtered  between  0.1  to  0.4  Hz. 
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Figure  38.  KKAR  receiver-function  record  sections.  The  receiver  functions  are  aligned  on 
the  minimum-phase  time  and  have  been  filtered  between  0.1  to  0.4  Hz. 
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The  moveout  of  the  radial  component  reverberations  and  the  existence  of  coherent  energy  on 
the  tangential  components  receiver  functions  seem  to  imply  that  the  earth  structure  below  the 
arrays  is  not  laterally  homogeneous,  as  we  expected  from  the  array  measurements.  In  order  to 
understand  the  earth  structure  we  utilize  a  simple  forward  modeling  approach.  We  start  with  a 
simple  flat-layer  velocity  model  determined  from  inverting  a  cumulative  stack  of  the  radial 
receiver  functions.  Using  this  model,  we  systematically  vary  the  strike  and  dip  of  the  Moho  to 
compute  a  suite  of  synthetic  receiver  functions.  The  synthetic  radial  and  tangential  record  sections 
are  then  quantitatively  compared  to  the  observed  record  sections.  In  this  manner  we  attempt  to 
solve  for  the  strike  and  dip  that  minimizes  the  difference  between  the  observed  and  computed 
record  sections. 

Results  using  this  grid  search  approach  for  the  orientation  of  the  converting  interface  below  the 
MKAR  array  were  inconclusive.  Our  analysis  revealed  that  the  Moho  below  the  MKAR  array  is 
47±6  km  deep.  For  the  KKAR  array,  the  analysis  suggested  a  Moho  at  a  depth  38±5  km  with  a 
strike  of  N30°E±20  and  a  dip  of  6°±10.  Results  from  these  exercises  indicate  that  the  near-array 
structure  is  quite  complicated  at  both  arrays,  and  more  research  is  required  into  methods  that  can 
account  for  near-array  effects  in  a  simple  and  straightforward  manner. 


7.2  Analysis  of  the  PcP  Phase 

In  another  attempt  to  gain  further  insight  into  near-array  structure,  we  examined  the  PcP  phase. 
At  far-regional  distances,  PcP  arrives  with  a  nearly  vertical  incidence  angle  and  any  moveout 
across  the  array  should  be  directly  related  to  the  near-array  structure.  Our  plan  was  to  incorporate 
results  from  the  PcP  analysis  with  the  receiver  function  results  to  provide  a  clearer  picture  of  the 
structure  below  the  array.  This  could  lead  to  a  better  understanding  of  the  heterogeneity  that  is 
important  at  far-regional  distances. 

Figure  38  shows  the  expected  slowness  and  travel-time  of  the  PcP  phase  for  far-regional 
distances.  Between  10  -  30°  distance,  PcP  arrives  with  phase  velocities  between  100  and  40  km/s, 
which  at  the  shorter  offsets  is  essentially  vertically  incident.  It  is  this  feature  that  we  wish  to  exploit 
to  help  better  understand  non-planar  structures  situated  below  the  arrays.  However,  identifying 
PcP  can  be  difficult  at  far-regional  distances  for  several  reasons,  including  weak  amplitudes  that 
are  due  to  the  impedance  at  the  core-mantle  boundary  and  potential  interference  from  S-wave 
arrivals  at  a  distance  of  ~20°. 

We  identified  56  events  from  our  far-regional  database  that  had  potentially  observable  PcP 
arrivals.  These  events  were  observed  between  12  -  26°  epicentral  distances  and  had  a  good  overall 
signal-to-noise  quality.  In  Figure  39  we  show  a  spectrogram  of  signal  windowed  around  the 
AK135-predicted  PcP  arrival  for  an  event  recorded  by  the  KKAR  array  at  a  distance  of  16.4°.  In 
this  particular  example  several  packets  of  energy  arrive  within  a  20-second  window  of  the 
predicted  PcP  arrival.  To  identify  these  signals,  we  applied  array-processing  methods  to  find 
arrivals  occurring  at  near  vertical  incidence. 
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Figure  39.  Behavior  of  the  PcP  phase  at  far-regional  distances.  Left:  Slowness  computed  for 
PcP  for  a  source  at  the  Earth’s  surface  using  the  AK135  reference  model;  phase  velocity 
at  10°  (~100  km/s)  and  30°  (40  km/s)  are  labeled.  Right:  The  travel-time  of  PcP  for  source 
depths  of  0,  10,  20,  40  and  80  km. 

Figure  40  shows  the  slowness  vespagrams  for  the  event  in  Figure  39.  In  this  case,  we  used 
phase-coherence  semblance  stacking  (described  in  a  previous  section)  to  measure  the  phase 
velocities.  Most  of  the  energy  in  the  window  arrives  at  high  phase  velocities  indicating  a  steep 
incidence  to  the  array.  The  predominant  arrival  occurs  at  ~13  seconds  with  a  velocity  of  ~45  km/s. 
This  is  preceded  6  seconds  earlier  by  a  weaker  arrival  with  a  similar  velocity  (both  are  labeled  in 
Figure  40).  Such  high  phase  velocities  limit  the  potential  phases  these  arrivals  represent,  ruling  out 
S  arrivals  or  surface  waves. 

We  processed  all  56  events  in  our  data  set  for  PcP  arrivals.  In  order  to  confidently  identify 
these  phases,  we  required  observations  from  several  events  over  a  range  of  distances  that  exhibited 
high  phase  velocities.  Out  of  the  56  events,  we  were  able  to  find  9  with  potential  PcP  arrivals.  In 
some  cases,  possible  PcP  precursors  were  also  observed,  and  slower  arriving  energy  in  the  same 
time  window  was  interpreted  as  potential  PcP  converted  phases  from  structure  near  the  array. 
However,  further  analysis  would  be  needed  to  verify  these  hypotheses. 
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Figure  40.  Spectrogram  of  a  seismogram  window  around  the  predicted  PcP  phase  (zero 
time).  The  seismic  event  was  recorded  by  the  KKAR  array  at  a  distance  of  16.4°  and  the 
central  array  element  is  shown.  Red  circle  highlight  potential  PcP  arrival. 
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Figure  41.  Phase-weighted  semblance  stacking  for  PcP  arrivals.  Time  axis  is  the  same  as 
Figure  39  .  Data  were  filtered  between  0.8  and  1.75  Hz  prior  to  processing.  Potential  PcP 
arrivals  with  high  phase  velocities  are  labeled. 
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In  order  to  attempt  to  deduce  information  on  potentially  dipping  structure  below  the  arrays  we 
cross-correlated  the  PcP  arrivals  across  the  individual  arrays  to  determine  the  time  delay  of  the 
phase.  Since  PcP  is  near  vertical  at  short-offsets,  it  was  assumed  that  the  pattern  in  the  time  delays 
would  be  directly  related  to  orientation  of  the  non-planar  structures  below  the  array.  However,  for 
7  of  the  9  events  the  cross-correlation  delay  times  were  0.0  for  most  array  elements,  suggesting 
that  either  the  structure  below  the  arrays  is  planar  or  our  assumption  concerning  PcP  was  incorrect. 
For  the  other  two  events,  the  cross-correlation  delay  times  could  not  be  fit  with  a  planar  surface. 

In  summary,  the  near-array  structure  below  MKAR  and  KKAR  appears  to  be  quite 
complicated.  It  has  a  strong  affect  on  back-azimuth  residuals  (±20°)  and  receiver  functions  show 
large  amounts  of  tangential  energy.  Our  attempts  to  quantify  the  structure  using  a  single  dipping 
interface  were  not  completely  successful.  At  MKAR  we  were  not  able  to  obtain  a  reliable  solution 
with  such  a  model  from  either  the  polarization  analysis,  receiver  function  analysis,  or  by 
examining  PcP  arrivals.  At  KKAR  the  polarization  and  receiver- function  analysis  suggest  a  Moho 
that  dips  to  the  northeast  at  6  -  10°. 
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Appendix  1:  Altering  Small- Aperture  Array  Geometry  to  Improve  Phase-Velocity 
Estimation 

A  major  focus  of  this  project  was  the  development  of  analysis  methods  to  process  far-regional 
arrivals  observed  on  small-aperture  arrays.  By  analyzing  the  array  response  functions,  we  know 
that  the  small  array  apertures  of  MKAR  and  KKAR  produce  a  lack  of  precision  in  slowness  and 
back-azimuth  estimates,  and  the  inter-sensor  spacing  creates  wave-number  aliasing  that  reduces 
accuracy  in  the  slowness  estimate.  In  this  appendix  we  show  that  significant  improvements  in 
slowness  resolution  can  be  gained  by  analyzing  data  from  a  sub-array  from  the  current  MKAR  and 
KKAR  configuration,  or  by  adding  a  single,  optimally  placed  array  element  at  some  distance  away 
from  the  center  of  the  array. 

For  example,  Figure  A-l  shows  the  current  configuration  of  the  KKAR  array  (black  triangles), 
which  has  nine  elements  and  a  total  aperture  of  ~3.5  km.  At  far-regional  distances,  the  P  arrivals 
travel  across  the  array  at  phase  velocities  between  8.0  -  12.0  km/s,  requiring  only  0.27  -  0.44 
seconds  to  traverse  the  full  aperture  of  the  array,  and  the  travel  time  between  adjacent  array 
elements  is  even  smaller.  While  the  small  differential  travel  time  across  the  array  allows  for  high 
signal  coherence  between  array  elements,  it  also  negatively  affects  the  quality  of  beam- forming 
results,  which  involve  time  shifting  and  stacking  to  find  the  directional  parameters.  This  makes 
distinguishing  a  Pnl  arrival,  or  a  depth  phase,  from  a  P410  arrival,  extremely  difficult  with 
slowness  (or  phase  velocity)  measurement  alone.  The  lack  of  precision  and  resolution  is  a 
consequence  of  the  small  time  shifts  required  between  array  elements  to  span  the  phase  velocity 
range  of  8.0  -  12.0  km/s  during  the  beam-forming  process.  The  small  time  shifts  do  not 
significantly  affect  the  power  or  coherency  in  the  beam  stack,  and  hence  the  confidence  regions  in 
arrival  measurements  span  a  large  range  of  slowness  and  back-azimuths. 
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Figure  A-l.  Black  triangles  depict  the  current  configuration  of  the  KKAR  array;  the  red 
triangles  are  hypothetical  stations  positioned  at  4,  8  and  12  km  from  the  center  of  the 
current  array. 

To  demonstrate  the  usefulness  of  altering  small-aperture  array  geometry,  we  performed 
vespagram  analysis  of  synthetic  P-wave  seismograms  generated  using  various  array 
configurations.  We  modeled  the  wavefield  for  a  point-source  dislocation  in  Iran  as  it  would  be 
observed  at  KKAR  at  a  distance  of  16°  and  a  back  azimuth  of  225°.  We  examined  the  vespagrams 
results  for  several  cases,  including  the  full  array,  and  sub-arrays  with  or  without  additional  outer 
elements.  This  hypothetical  example  represents  a  best-case  scenario,  since  the  synthetic  data  do 
not  contain  noise,  and  the  P  arrival  is  a  single  well-defined  pulse  (Figure  A-2). 

Array  Configuration  1:  Vespagram  Analysis  using  Sub-arrays 

In  Figure  A-3  we  show  the  results  of  applying  4th-root  vespagrams  analysis  to  two 
configurations  of  the  current  KKAR  array.  Figure  A-3a  shows  the  results  using  all  nine  elements 
of  the  KKAR  array,  and  these  can  be  compared  to  the  results  in  Figure  A-3b,  which  were 
computed  using  only  the  five  outer  array  elements  of  the  KKAR  array.  For  this  synthetic  example, 
the  vespagram  in  Figure  A-3b  is  more  peaked.  This  occurs  because  the  beam  power  built  up  by  the 
contribution  of  the  inner  elements  to  the  waveform  stack  changes  little  for  different  phase-velocity 
time-shifts,  due  to  spatial  closeness  of  the  array  elements.  In  essence,  using  all  the  inner  elements 
over-weights  the  stack  and  produces  the  phase-velocity  smearing  seen  in  Figure  A-3a.  We  note 
that  both  vespagram  analyses  produced  the  correct  phase  velocity  of  8.5  km/s  for  the  initial  P 
arrival.  The  sharpness  of  the  peaks  in  phase  velocity  is  best  illustrated  with  slices  taken  through  t  = 
0.825  s  in  each  vespagram;  these  are  shown  in  Figure  A-5  for  all  the  examined  cases. 
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Figure  A-2.  Synthetic  array  gather  for  the  current  configuration  of  KKAR  and  three 
additional  hypothetical  stations  (labeled  as  ??01,  ??02,  and  ??03).  Waveforms  were 
generated  from  a  point  source  in  Iran,  at  a  distance  of  16°  and  back  azimuth  of  225°  from 
KKAR.  The  waveforms  are  aligned  based  on  their  moveout  with  respect  to  the  closest 
station  to  the  event  (KK07  in  this  case).  The  P  arrival  and  the  vespa  analysis  window  are 
indicated  near  the  top  of  the  plot. 


Figure  A-3.  Vespagram  analysis  of  the  P  arrivals  in  the  synthetic  waveforms  in  Figure  A-2. 
a)  results  computed  using  all  nine  elements  of  the  array  (shown  to  the  right);  b)  results 
computed  using  only  the  outermost  five  elements  of  the  array;  the  phase  velocity  is  more 
sharply  defined. 
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Array  Configuration  2:  Vespagram  Analysis  with  Additional  Elements 

For  this  case  we  illustrate  the  effect  that  increasing  the  overall  array  aperture  has  on  the 
performance  of  the  vespagram  method.  By  increasing  the  array  aperture,  the  travel  time  across  the 
array  for  the  body-wave  arrivals  also  increases,  which  subsequently  improves  the  phase  velocity 
estimate.  Adding  a  station  outside  the  current  array  at  carefully  chosen  azimuths  can  also 
significantly  improve  phase-velocity  resolution  for  events  in  regions  of  interest.  For  example,  in 
Figure  A-4  we  show  the  vespagrams  computed  with  the  outer  ring  of  elements  at  KKAR  and  a 
single  array  element  at  4,  8,  or  12  km  from  the  center  of  the  array  along  the  great-circle  path  from 
the  event  to  the  southwest  (red  stations  in  Figure  A-l).  The  phase-velocity  peaks  become  sharper 
as  the  distance  of  the  additional  element  increases  from  the  center  of  the  array.  Adding  a  station  4 
km  from  the  array  center  (Figure  A-4a)  removes  a  significant  portion  of  the  velocity  smearing  seen 
in  the  original  analysis  in  Figure  A-3a.  The  least  amount  of  phase-velocity  smearing  is  seen  when 
a  station  is  added  at  8  (Figure  A-4b)  or  12  km  (Figure  A-4c)  from  the  array  center.  However,  in  a 
real  world  scenario,  site  noise  and  structural  heterogeneity  around  the  array  would  limit  the 
maximum  beneficial  distance  of  the  added  station,  since  maintaining  signal  coherence  across  the 
array  is  also  very  important  for  array  performance.  We  also  note  that  expanding  the  array  aperture 
without  adding  other  interior  array  elements  can  cause  unwanted  spatial  aliasing.  However,  both  of 
the  scenarios  included  here  produce  a  vespagram  with  little  smearing  at  high  phase  velocities. 
Figure  A-5  shows  the  narrow  (well-resolved)  peaks  in  phase  velocity  produced  by  the  adding 
single  elements  at  various  distances. 
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Figure  A-4.  Vespagrams  for  the  P  arrival  in  Figure  A-2,  using  the  five  outer  elements  of  the 
array  and  an  additional  hypothetical  station  at  a)  4  km;  b)  8  km;  and  c)  12  km  from  the 
center  of  the  array  and  along  the  great-circle  path  from  the  event.  As  the  distance 
between  the  station  and  array  increases,  the  P  phase  velocity  becomes  significantly  easier 
to  pick  with  confidence. 
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Figure  A-5.  Profile  comparisons  of  the  vespagrams  shown  in  Figures  A-3  and  A-4.  The 
color-coded  lines  are  profiles  taken  through  the  t  =  0.825  s  in  each  vespagram.  While  all 
profiles  have  their  absolute  peak  at  the  correct  phase  velocity  (8.5253  km/s),  the 
phase-velocity  resolution  is  improved  by  increasing  the  array  aperture. 

In  summary,  the  small  aperture  of  regional  arrays  such  as  KKAR  limits  the  velocity  resolution 
for  events  observed  at  distances  between  13  -  30°.  At  these  distances,  multiple  arrivals  occur  in 
short  time  windows,  and  the  velocity  resolution  becomes  critically  important  to  differentiate 
between  them.  While  our  analysis  was  limited  to  synthetic  data,  it  appears  that  the  resolution  in 
phase  velocity  can  be  improved  by  simply  limiting  the  number  of  inner  array  elements  used  in 
vespagram  analysis.  However,  the  greatest  improvement  occurs  when  stations  are  added  outside  of 
the  current  array  configuration.  While  we  did  not  determine  the  optimal  array  configuration  to 
analyze  all  intermediate-distance  events,  the  examples  shown  here  suggest  that  adding  3  to  5  array 
elements  at  8  -  12  km  distance  from  the  array  center  would  significantly  improve  array  analysis 
methods  for  events  at  this  distance  range. 
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